### KARLSRUHER INSTITUT FÜR TECHNOLOGIE (KIT) SCHRIFTENREIHE DES INSTITUTS FÜR TECHNISCHE MECHANIK

BAND 36

JENS BURGERT

On direct and inverse problems related to longitudinal impact of non-uniform elastic rods

36

Jens Burgert

## On direct and inverse problems related to longitudinal impact of non-uniform elastic rods

Jens Burgert

**On direct and inverse problems related to longitudinal impact of non-uniform elastic rods**

### **Karlsruher Institut für Technologie Schriftenreihe des Instituts für Technische Mechanik**

Band 36

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

## **On direct and inverse problems related to longitudinal impact of non-uniform elastic rods**

by Jens Burgert

Dissertation, Karlsruher Institut für Technologie KIT-Fakultät für Maschinenbau

Tag der mündlichen Prüfung: 07. Dezember 2020 Hauptreferent: Prof. Dr.-Ing. Wolfgang Seemann Korreferent: Prof. Dr. Dr. h.c. Peter Hagedorn

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2021 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 1614-3914 ISBN 978-3-7315-1087-1 DOI 10.5445/KSP/1000129237

# **Danksagung**

Die vorliegende Arbeit entstand während meiner Tätigkeit als wissenschaftlicher Mitarbeiter am Institut für Technische Mechanik (ITM), Abteilung Dynamik/Mechatronik des Karlsruher Instituts für Technologie (KIT).

An erster Stelle gilt mein Dank Herrn Prof. Dr.-Ing. Wolfgang Seemann für die fachliche Betreuung. Ich danke ihm für die vertrauensvolle Zusammenarbeit, die wissenschaftliche Freiheit und die stete Unterstützung und Förderung. Hervorheben möchte ich seine Unterstützung beim Aufbau internationaler Kooperationen mit anderen Universitäten, von denen ich fachlich und persönlich sehr profitiert habe.

Weiterhin bedanke ich mich bei Prof. Dr. Dr. h.c. Peter Hagedorn für die Übernahme des Korreferats und den damit verbundenen Anregungen zu meiner Arbeit. Herrn Prof. Dr. mont. Christoph Kirchlechner danke ich für die Übernahme des Prüfungsvorsitzes.

Ebenfalls danke ich Herrn Prof. Dr.-Ing. habil. Alexander Fidlin für die stets ehrlichen Fragen und Anmerkungen in den Seminaren. Herrn Prof. Dr.-Ing. Carsten Proppe danke ich für die angenehme Zusammenarbeit am Institut.

Den emeritierten Professoren Prof. Dr.-Ing. Dr. h.c. Jörg Wauer, Prof. Dr.-Ing. Walter Wedig und Prof. Dr.-Ing. Dr. h.c. Jens Wittenburg danke ich für viele interessante Gespräche. Ihre Teilnahme am Institutsleben ist eine große Bereicherung. Ich hoffe, dass Sie dem wissenschaftlichen Nachwuchs noch viele Jahre zur Seite stehen.

Die Forschungsaufenthalte an der Universität Tampere in Finnland und der Universität Uppsala in Schweden haben der vorliegenden Arbeit entscheidende Impulse gegeben. Großen Anteil am Gelingen der Aufenthalte haben die Professoren Prof. Dr. Erno Keskinen und Prof. Dr. Bengt Lundberg. Ihnen danke ich von Herzen.

Des Weiteren bedanke mich bei allen Kolleginnen und Kollegen für die wundervollen Jahre am Institut. Neben den fachlichen Diskussionen werden mir vor allem das kollegiale Miteinander, die Hilfsbereitschaft und die persönlichen Gespräche in positiver Erinnerung bleiben. Hervorheben möchte ich meinen Zimmerkollegen Jimmy Aramendiz, mit dem ich über jedes Thema philosophieren konnte und mit dem ich auch abseits der Tätigkeit am Institut viele schöne Stunden verbringen durfte.

Ebenfalls danke ich allen Abschlussarbeitern und Hilfswissenschaftlern, deren Betreuung mir sehr viel Freude bereitet hat und von deren Arbeit und Fragen ich sehr profitiert habe.

Ganz besonders möchte ich mich bei meiner Familie bedanken. Ihr konstanter Rückhalt und ihre Unterstützung in jeglicher Hinsicht haben mich tief beeindruckt und berührt. Durch sie habe ich nie aus den Augen verloren, was die wirklich wichtigen Dinge im Leben sind. Dabei möchte ich besonders meinen Eltern für ihre Hingabe und ihre bedingungslose Liebe danken.

Zu guter Letzt gilt mein größter Dank meiner lieben Frau Kristin, die mich stets motiviert und unterstützt hat, insbesondere in dem turbulenten letzten Jahr. Ich danke ihr vor allem dafür, dass sie mein Leben auf eine so wundervolle Art und Weise bereichert.

Karlsruhe, den 10. Dezember 2020 Jens Burgert

# **Zusammenfassung**

Das Ziel der Arbeit ist die Entwicklung einer Methode zur Untersuchung inverser Probleme in aufeinandertreffenden Stäben. Dabei soll die Impedanz des auftreffenden Stabs ermittelt werden, sodass eine vorgegebene Stoßkraft erzeugt wird (*Impedanzproblem*). Die Methode ist für viele experimentelle und praktische Anwendungen relevant. Beim Schlagbohrverfahren stellt sich beispielsweise die Frage, wie die Impedanz des auftreffenden Stabs gewählt werden kann, sodass die einlaufende Kraftwelle bezüglich der Effizienz optimiert wird.

Das *Impedanzproblem* wird nur von wenigen Forschern behandelt. Der vielversprechendste Ansatz löst das inverse Problem mithilfe einer Integralgleichung. Dabei werden mehrere Näherungsmethoden angewandt, um die Impedanz des auftreffenden Stabs zu ermitteln. Im Rahmen dieser Arbeit wird eine Methode entwickelt, die auf der Idee differentieller Methoden basiert. Die Methode nutzt die Struktur des Wellenausbreitungsproblems aus. Hierfür werden die Stäbe durch Elemente mit stückweiser konstanter Impedanz diskretisiert. Zunächst wird die Impedanz des ersten Elements ermittelt, welches direkt auf die gestoßene Stange auftrifft. Alle weiteren Elemente beeinflussen die Stoßkraft zeitverzögert. Je weiter die Elemente von der Stoßfläche entfernt sind, desto später ist ihr Einfluss sichtbar. Zum Zeitpunkt des erstmaligen Einflusses des jeweiligen Elements wird die zugehörige Impedanz so gewählt, dass die vorgegebene Stoßkraft erreicht wird. Dadurch werden die diskretisierten Impedanzen des auftreffenden Stabs rekursiv bestimmt.

Für stückweise konstante Impedanzfunktionen liefert die entwickelte Methode exakte Lösungen in geschlossener Form. Des Weiteren wird eine Bedingung hergeleitet, die angibt, ob eine Lösung zu dem inversen Problem existiert. Es wird gezeigt, dass jedes zusätzliche Element, dessen Impedanz sich von der Impedanz des vorigen Elements unterscheidet, die Wahrscheinlichkeit der Existenz einer Lösung verringert.

Das zugrundeliegende 1D-Modell wird experimentell validiert, indem die Simulationsergebnisse mit Messungen verglichen werden. Es wird insbesondere gezeigt, dass das Modell auch für komplexe Geometrievariationen der Stäbe gültig ist.

# **Abstract**

Impacting rods are used in many practical and experimental applications. The analysis problem of finding the time-dependent impact force for given properties of the impacting rods and known initial velocities is thoroughly studied in literature. However, very little literature has been published which addresses its inverse synthesis problem: Find the location-dependent impedance function of the impacting rod that generates a prescribed impact force.

In this contribution, a method that solves the synthesis problem is developed. The impedances of the linear elastic rods are discretized by elements of piecewise constant impedance. This enables making use of the analytical solution which exists for these elements. The method utilizes that each element of the unknown impactor influences the prescribed force at different time instants for the first time. At this specific time, the impedance of the corresponding element is adjusted so that the actual force matches the prescribed force. This procedure is repeated recursively until the impedances of all impactor elements are determined. For piecewise constant impedance functions, the developed method delivers exact solutions in closed-form. Moreover, a condition is derived which states if a physically meaningful solution exists for the considered element. Thus, it is known if and how long a solution of the inverse problem exists. The working principle of the method is demonstrated with several examples, including the application in percussion drilling.

Finally, the underlying 1D model is validated experimentally. To this end, a single-hit test rig is set up. The comparison of simulation versus experimental results yields good accordance both for standard and complex geometries of the impacting rods.

# **Contents**



# **1 Introduction**

Wave propagation phenomena are present in anyone's daily routine. Every time we talk to someone, pressure waves are generated which propagate through the air to the ears of the listener. The sound is created by inhaled air which is sent back out of the mouth. On the way from the lung to the mouth, the airflow is shaped by the vocal cords (excitation) and the cross sectional area variation of the vocal tract. Humans intuitively choose the appropriate excitation and cross section area of the vocal tract while talking. However, it is of interest to identify which vocal tract cross section area leads to the corresponding sound. The direct problem of determining the pressure at the lips for given excitation and cross section area of the vocal tract is very well understood. As opposed to this, the inverse problem of finding the vocal tract's cross section area for given excitation and pressure measurement at the lips is much more demanding.

Typical questions that arise when dealing with inverse problems are:


Another example of an inverse wave propagation problem is the determination of the characteristic rod impedance function by applying a force impulse at one end of the rod and measuring the impulse response. Similar questions are found in transmission line analysis and geophysics. Mathematically, the presented scattering problems are described by the same hyperbolic system of partial differential equations.

The inverse scattering problem of determining the rod's impedance is of little practical interest as the impedance is either known or determined by other means. However, practical applications emerge if the system is extended by an impacting rod. For example, when using percussion drilling to generate blast holes for mining, an actuated piston

transfers its kinetic energy to the rod during the impact. The induced force waves propagate along the rod until they reach the ground. The efficiency of the drilling process largely depends on the shape of the incident force waves at the ground. This motivates the development of a method that allows to prescribe the force over time at any axial rod position and to determine the corresponding impedance function of the impacting rod. A new method is necessary since the classical approaches of the scattering problem are not directly applicable to the impacting rod problem.

### **1.1 State of Research**

In Fig. 1.1, an overview of the state of research is given. The classical 1D scattering problem consists of applying an impulse 𝐼(0, 𝑡) at the boundary of the system and measuring the impulse response 𝑂(0, 𝑡). As the input and output functions are known, the properties of the system may be reconstructed. The wave supporting system is described by two partial differential equations with coefficients <sup>𝛾</sup>(𝑥) and <sup>𝜁</sup>(𝑥) [46]. Physically, the coefficients refer to the distributed parameters of the system. Depending on the interpretation of the coefficients, the model is applicable in many fields of engineering. Some examples are the determination of the earth properties in geophysics [8, 10, 63, 64], the determination of the vocal tract shape by acoustical measurements at the lips [4, 65, 68 – 70] and the analysis of non-uniform transmission lines [12, 15, 51, 52].

There are different names for the partial differential equations in literature [11] depending on the parametrization. The *telegrapher's equations* result from the original equations by carrying out a coordinate transformation and introducing the impedance. After normalizing the *telegrapher's equations*, the equations in the form of *Schrödinger equations* arise by introducing two potentials. For expressing the variables 𝐼(𝑥, 𝑡) and 𝑂(𝑥, 𝑡) in terms of left and right propagating waves, no established term exists in literature. However, it is sometimes called *two component wave system*. Since the equations may be transformed into the original system, they all rely on the same mathematical foundation.

As many engineering fields are involved in finding solutions to the inverse scattering problem, the list of published literature is large. A comprehensive overview of many scattering problems is presented in [30] and [11]. The idea of solving the scattering problem using perturbation theory [54] has not been developed further since it has several disadvantages. The most serious ones are that the method does not make use of the structure of the problem and that the boundary conditions have to be known. The most promising methods are either assigned to differential methods or integral


**Figure 1.1:** Overview of the state of research and the content of the thesis.

equations. Historically, the first complete solution of the scattering problem has been documented by Gelfand and Levitan in terms of an integral equation. Subsequently, many similar integral equations have been published, e.g. in references [3, 68, 74]. They mainly differ in the input/output pair that is used in the formulation. In general, the computational costs of the integral equations is 𝒪(𝑛¯ 3 ), with 𝑛¯ being the number of discretized time intervals. However, by making use of the special structure of the problem, the costs are reduced to 𝒪(𝑛¯ 2 ) [5].

The differential methods directly exploit the knowledge of the wave propagation process. As the material properties of the medium only vary in one direction, the medium is interpreted as layers of constant material properties. The further the layers are away from the force excitation the longer it takes to reach them. This also applies to the measured signal which is influenced by the layers at different moments for the first time. These observations are utilized in differential methods which determine the physical properties of the medium layer for layer [14, 62, 69, 71]. For a continuous change of the layers, the method is usually called *layer peeling method* [11]. The method with discrete changes of the layers is known as the *dynamic deconvolution method* [60] or as the *downward continuation method* [14].

The presented methods are also applicable to the longitudinal wave propagation in elastic rods. The classic scattering problem of determining the impedance of the rod by its impulse response is of little practical interest because the impedance is either known or determined by other means. However, the methods for the scattering problem may be extended to impacting rods where the impedance of one rod is determined so that a desired impact force is generated (*impactor synthesis problem*). This problem is motivated by theoretical work of Lundberg and Collet who have shown that the efficiency of the percussion drilling process is increased if the incident wave is of exponential shape [49]. Lundberg and Lesser have extended the integral equations of Gopinath and Sondhi [68] to solve the *impactor synthesis problem* in terms of an integral equation [47]. Dutta has suggested to calculate the rod's impedance by systematically solving the direct problem and determining the solution which fits best to the desired impact force [24].

The *impactor synthesis problem* also occurs in dynamic material testing using a Split Hopkinson Pressure Bar apparatus. The test rig was invented by Kolsky [40] and has been applied to determine the material properties of a test specimen at high strain rates [73]. The setup consists of three bars (striker, incident, transmitter) and a test specimen which is located between the incident and the transmitter bar. After the striker hits the incident bar, stress waves propagate along the incident bar toward the test specimen. At the test specimen, the stress waves are partly reflected and transmitted. The incident, reflected, and transmitted waves are recorded by strain gauges that are attached to the incident and transmitter bar. For the evaluation of the detected waves, several methods have been proposed, e.g. in references [29, 33, 55, 57].

It might be attempted to generate an incident wave that produces a constant strain rate in the specimen. An evident difficulty is that this wave depends on the unknown properties of the specimen. One possibility to fulfill this requirement is to optimize the geometry of the striker [22, 26, 45, 75] accordingly. Another approach is to introduce a shaping bar between striker and incident bar [28]. So far, the optimized striker geometries are designed by experience and knowledge of the solution of the direct problem. However, the determination of the geometry by solving the corresponding inverse problem promises to deliver faster and more accurate results.

Apart from the two presented approaches, no other methods are known in literature which solve the *impactor synthesis problem*. Specifically, no methods are known which are based on the idea of differential methods to recursively determine the rod's impedance layer by layer.

The quality of the approaches strongly relies on the validity of the 1D model. Pochhammer [58] and Chree [21] independently formulated the exact solution of the three-dimensional wave propagation problem in an infinitely long cylinder. The solution relates the phase velocity of the longitudinal waves to the wavelength of its propagating modes (dispersion relation) [6]. The dispersive effect increases the greater the ratio between rod diameter and wavelength is [27]. To correct these effects, several methods have been established, e.g. in references [25, 53].

Since the dispersive effects are not represented by the 1D model, it has to be examined under which circumstances the 1D model delivers accurate results. It has been proved that the 1D results yield good accordance with their corresponding experimental results, e.g. in reference [56]. However, only standard geometries (constant, conical, etc.) have been studied. It has not been systematically examined yet if a complex variation of the rod geometry, including jumps in cross section area, delivers good accordance with the experimental results.

## **1.2 Thesis Purpose and Structure**

The literature research reveals that there are very few methods which cover the *impactor synthesis problem*. One approach is to systematically optimize the rod's impedance by solving the direct problem. This approach has several drawbacks: On the one hand, the method does not guarantee that the solution converges to the exact solution and on the other hand, the computational costs are very high. A more promising approach has been developed by Lundberg and Lesser [47] based on the integral equations of Gopinath and Sondhi [68]. This method provides a solution in two steps: Initially, a deconvolution is carried out to determine the impulse response of the unknown rod. Subsequently, the impulse response is substituted into an integral equation which is finally evaluated for the determination of the unknown rod impedance. However, several approximation methods are applied to obtain the rod's impedance. Firstly, the impulse response of the given rod, which is part of the convolution equation, is not always known analytically and has to be approximated. Secondly, approximation methods are applied to deconvolve the convolution equation. And thirdly, the integral equation is solved approximately.

As compared to the evaluation of integral equations, a differential-based method, that exploits the structure of the problem, merely has to approximate the unknown rod's impedance by piecewise constant layers. However, the existing methods to solve the scattering problem are not expandable to the *impactor synthesis problem* straightforwardly.

The **main purpose** of the thesis is to develop a method that solves the *impactor synthesis problem* based on the ideas of the differential methods of the scattering problem. A **further aim** is to validate the underlying 1D model experimentally for an impact with complex rod geometries.

In the thesis, the setups presented in Fig. 1.2 are investigated. In Setup 1, a rod with constant cross section area is impacting a rod with exponentially increasing cross section area whose right end is fixed. Setup 2 is a textbook example where a rod consisting of two cross section area sections is hitting a rod with a constant cross section area. Setup 3 is more complex, including variations in geometry and the consideration of external friction forces. Setup 4/5 only differ in the prescribed impact force <sup>𝑁</sup><sup>I</sup> that shall be generated by the impacting rod whose impedances are determined accordingly. In Setup 6, a rod with two steps in cross section area is hitting a rod with a constant cross section area. Setup 7 is similar to Setup 6 but with more complex cross section area

**Figure 1.2:** Investigated setups.



**Figure 1.3:** Overview of the content of the chapters.

variations of the impacting rod. The setups are investigated by different means. They are either investigated by applying the modal expansion method, the wave propagation method, the recursive method or experimentally (see Fig. 1.3). The main purpose of the thesis is to develop the recursive method which solves the *impactor synthesis problem*. The recursive method is based on the wave propagation method. The modal expansion method verifies the wave propagation method whereas the experimental results validate the underlying 1D model. This motivates the structure of the following chapters:

**Chapter 2:** In this chapter, the fundamentals that are used in the following chapters are presented.

**Chapter 3:** The wave propagation method for the calculation of 1D wave propagation problems in non-uniform rods is introduced. The range of applications is demonstrated by applying the method to Setup 1 and Setup 2. The method is verified by a comparison of results obtained with the modal expansion method versus results obtained with the wave propagation method applied to Setup 3.

**Chapter 4:** The fourth chapter covers the new recursive method. The presentation of the main idea is followed by its mathematical description. A condition is derived to assess the existence of a solution of the inverse problem for every step. Afterward, the working principle of the recursive method is presented by employing Setup 4/5. Setup 4 only differs from Setup 5 in the prescribed impact force <sup>𝑁</sup><sup>I</sup> . Depending on the impact force, a solution of the inverse problem exists (Setup 4) or not (Setup 5). Finally, the method is applied in percussion drilling.

**Chapter 5:** The experimental setup is presented. Several experiments are conducted using Setup 6 to determine the influence factors on the measurements and to examine the repeatability of the results. Subsequently, the simulation and experimental results obtained with Setup 6 and Setup 7 are analyzed and compared. Thus, the underlying 1D model is validated.

**Chapter 6:** In the last chapter, the thesis is concluded.

# **2 Fundamentals**

In this chapter, the fundamentals, which are applied in the subsequent chapters, are presented. Initially, the governing partial differential equation (PDE) for the 1D longitudinal wave propagation in rods is derived. It is shown that the PDE can be solved analytically for a special class of cross sectional area variations by applying the modal expansion method. However, the respective results are not expressable in closed-form but with an infinite series of eigenfunctions and eigenfrequencies. The modal expansion method reveals two drawbacks: Firstly, the analytical solution applies only to very few cross sectional area variations. And secondly, many series terms have to be considered to approach the exact solution. The solutions obtained by the modal expansion method are used to verify the results of the wave propagation method whose main idea is to subdivide the impacting rods into a large number of elements of piecewise constant impedance. For these elements, the wave equation can be solved easily. Its D'Alembert solution is discussed in detail. Especially, the transition and boundary conditions are considered as they are essential to solve the inverse problem.

### **2.1 Derivation of the Governing PDE**

In Fig. 2.1a, a schematic sketch of a rod is depicted. The mass density 𝜌, the Young's modulus 𝐸 and the cross section area 𝐴 only vary in 𝑥-direction. The axial displacement of the particles at position 𝑥 is described by 𝑢(𝑥, 𝑡). It is assumed that the centroid of the cross section area lies on a straight line perpendicular to the cross section area [36]. Applying Newton's second law on the free body diagram of a rod segment (see Fig. 2.1b) leads to

$$\frac{\partial N(\mathbf{x},t)}{\partial \mathbf{x}} = \rho(\mathbf{x})A(\mathbf{x})\frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} \tag{2.1}$$

**Figure 2.1:** (a) Schematic sketch of a rod and (b) the free body diagram of a rod segment for the derivation of the governing PDE.

where 𝑁(𝑥, 𝑡) is the normal force, positive in tension [35]. The particle velocity

$$w(\mathbf{x}, t) = \frac{\partial u(\mathbf{x}, t)}{\partial t} \tag{2.2}$$

is obtained by differentiating the displacement with respect to time whereas the axial strain

$$
\varepsilon(\mathbf{x}, t) = \frac{\partial u(\mathbf{x}, t)}{\partial \mathbf{x}}\tag{2.3}
$$

is obtained by differentiating the displacement with respect to space [50]. Inserting the force relation

$$N(\mathbf{x}, t) = A(\mathbf{x}) \sigma(\mathbf{x}, t) \tag{2.4}$$

together with the relation between axial stress <sup>𝜎</sup>(𝑥, 𝑡) and strain <sup>𝜀</sup>(𝑥, 𝑡) (Hooke's law)

$$
\sigma(\mathbf{x},t) = E(\mathbf{x})\varepsilon(\mathbf{x},t) = E(\mathbf{x})\frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}}\tag{2.5}
$$

into Eq. (2.1) yields the governing PDE

$$\frac{\partial}{\partial \mathbf{x}} \left[ E(\mathbf{x}) A(\mathbf{x}) \frac{\partial u(\mathbf{x}, t)}{\partial \mathbf{x}} \right] = \rho(\mathbf{x}) A(\mathbf{x}) \frac{\partial^2 u(\mathbf{x}, t)}{\partial t^2} \tag{2.6}$$

for the longitudinal wave propagation in rods. The model is equivalently describable by a hyperbolic system of linear partial differential equations

$$\begin{split} \frac{\partial N(\mathbf{x},t)}{\partial t} &= E(\mathbf{x})A(\mathbf{x}) \frac{\partial v(\mathbf{x},t)}{\partial \mathbf{x}}, \\ \frac{\partial v(\mathbf{x},t)}{\partial t} &= \frac{1}{\rho(\mathbf{x})A(\mathbf{x})} \frac{\partial N(\mathbf{x},t)}{\partial \mathbf{x}}. \end{split} \tag{2.7}$$

At this point, it is advantageous to introduce the impedance [32]

𝑍(𝑥) <sup>=</sup> 𝐴(𝑥)𝜌(𝑥)𝑐(𝑥) (2.8)

with

$$c(\mathbf{x}) = \sqrt{\frac{E(\mathbf{x})}{\rho(\mathbf{x})}}\tag{2.9}$$

being the wave propagation speed. The impedance of a rod is defined as the ratio between the force amplitude and its related velocity amplitude [36]. Thus, the impedance may be interpreted as the resistance of the elastic structure to the acting forces.

Introducing the coordinate

$$
\Theta(\mathbf{x}) = \int\_0^\mathbf{x} \frac{1}{c(\mathbf{x}')} \, \mathbf{d} \mathbf{x}' \tag{2.10}
$$

together with the notation 𝑣(𝑥(𝜃), 𝑡) <sup>=</sup> 𝑣¯(𝜃, 𝑡) etc. enables to transform the system (2.7) into [46]

$$\begin{split} \frac{\partial \bar{N}(\theta, t)}{\partial t} &= \bar{Z}(\theta) \frac{\partial \bar{v}(\theta, t)}{\partial \theta}, \\ \frac{\partial \bar{v}(\theta, t)}{\partial t} &= \frac{1}{\bar{Z}(\theta)} \frac{\partial \bar{N}(\theta, t)}{\partial \theta}. \end{split} \tag{2.11}$$

Therefore, the equations are describable solely in dependance of the characteristic impedance.

## **2.2 Analytical Solution of the Governing PDE**

In general, it is not possible to solve the governing PDE analytically. For a special class of cross sectional area variations and for homogeneous materials analytical solutions exist [1, 7, 43, 44, 59]. Guo and Yang [34] have shown that the cross sectional area variation has to be either square, exponential, square of hyperbolic cosine or square of cosine. The corresponding approach is named as the modal expansion method. Its main ideas are presented below.

For homogeneous materials, the governing PDE (Eq. (2.6)) is rewritten

$$\frac{\partial^2 u(\mathbf{x},t)}{\partial \mathbf{x}^2} + h(\mathbf{x}) \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}} - \frac{1}{c\_0^2} \frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} = 0, \text{ } h(\mathbf{x}) = \frac{\mathbf{d}A(\mathbf{x})/\mathbf{dx}}{A(\mathbf{x})}, \ c\_0^2 = \frac{E}{\rho}. \tag{2.12}$$

Inserting the exponential approach

$$u(\mathbf{x},t) = a(\mathbf{x})\mathbf{e}^{l(k\mathbf{x}-\omega t)}\tag{2.13}$$

with wave amplitude 𝑎(𝑥), wave number 𝑘 and angular frequency <sup>𝜔</sup> in Eq. (2.12) leads to two equations

$$\begin{aligned} \text{Im: } 2\frac{\text{da}(\mathbf{x})}{\text{d}\mathbf{x}} + h(\mathbf{x})a(\mathbf{x}) &= 0, \\ \text{d}^2 a(\mathbf{x}) &\quad \text{d}a(\mathbf{x}) &\quad \omega^2 \end{aligned} \tag{2.14}$$

$$\text{Re: } \frac{\text{d}^2 a(\mathbf{x})}{\text{d}\mathbf{x}^2} - k^2 a(\mathbf{x}) + h(\mathbf{x}) \frac{\text{d}a(\mathbf{x})}{\text{d}\mathbf{x}} + \frac{\omega^2}{c\_0^2} a(\mathbf{x}) = \mathbf{0} \tag{2.15}$$

for the imaginary (Im) and real (Re) part, respectively. From the first equation, 𝑎(𝑥) is determined

$$a(\mathbf{x}) = a\_0 \sqrt{\frac{A\_0}{A(\mathbf{x})}},\tag{2.16}$$

where <sup>𝑎</sup><sup>0</sup> and <sup>𝐴</sup><sup>0</sup> are the wave amplitude and cross section area at <sup>𝑥</sup> <sup>=</sup> 0. By eliminating 𝑎(𝑥) in Eq. (2.15), the equation

$$k^2 = \frac{\omega^2}{c\_0^2} - \underbrace{\left(\frac{1}{4}h(\mathbf{x})^2 + \frac{1}{2}\frac{\mathbf{d}h(\mathbf{x})}{\mathbf{d}\mathbf{x}}\right)}\_{b} \tag{2.17}$$

for 𝑘 is obtained. Analytical solutions are found if the wave number is independent of 𝑥 which means that 𝑏 has to be constant. This applies only to very few cross section areas, including the ones mentioned above. For those cross section areas the phase velocity

**Figure 2.2:** Phase speed for 𝜉 = 1 of the cross section areas that meet the requirement for an analytical solution (constant 𝑘).

$$c\_{\mathbb{P}} = \frac{\omega}{k} = \frac{c\_0}{\sqrt{1 - \frac{c\_0^2}{\omega^2}}} \tag{2.18}$$

depends generally on the angular frequency 𝜔 which can be seen in Fig. 2.2. Except from the curve related to quadratic cross sectional area variation, all the other curves are dispersive. Moreover, the exponential and square of hyperbolic cosine cross sectional area variations exhibit a cutoff frequency 𝜔<sup>c</sup> that is calculated with

$$
\omega\_c = c\_0 \sqrt{b}.\tag{2.19}
$$

Therefore, these cross section areas are mechanical high pass filters. In Tab. 2.1, the most important values related to the analytical solutions are summarized. These values are the cross section area variation 𝐴(𝑥), the constant 𝑏, the cutoff frequency <sup>𝜔</sup>c, the phase velocity <sup>𝑐</sup><sup>p</sup> and the wave amplitude <sup>𝑎</sup>(𝑥).


**Table 2.1:** Properties of the analytical solutions, based on [34].

## **2.3 Wave Equation**

The governing PDE (Eq. (2.6)) simplifies for homogeneous materials and constant cross section areas to the well known wave equation [23, 41]

$$\frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} = c^2 \frac{\partial^2 u(\mathbf{x},t)}{\partial \mathbf{x}^2}, \ c^2 = \frac{E}{\rho}. \tag{2.20}$$

The wave equation also applies to all physical values that are proportional to a partial derivative of 𝑢, especially, the velocity

$$w(\mathbf{x}, t) = \frac{\partial u(\mathbf{x}, t)}{\partial t},\tag{2.21}$$

the strain

$$
\varepsilon(\mathbf{x},t) = \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}},\tag{2.22}
$$

and the normal force

$$N(\mathbf{x}, t) = EA \frac{\partial u(\mathbf{x}, t)}{\partial \mathbf{x}}\tag{2.23}$$

as can be proven by partially differentiating the wave equation with respect to 𝑡 and 𝑥 [72].

### **2.3.1 D'Alembert Solution**

Initially, the wave propagation along an infinitely long rod without any transition or boundary conditions is considered. The two arbitrary functions 𝑓 (𝑥 <sup>−</sup> 𝑐𝑡) and 𝑔(𝑥 <sup>+</sup> 𝑐𝑡) with argument 𝑥 <sup>−</sup> 𝑐𝑡 and 𝑥 <sup>+</sup> 𝑐𝑡 solve the wave equation [2]. Due to the linearity of the PDE, the sum of the functions also solves the wave equation and the general solution of the displacement is given by

$$u(\mathbf{x},t) = f(\mathbf{x} - ct) + g(\mathbf{x} + ct). \tag{2.24}$$

The functions 𝑓 and 𝑔 describe waves traveling at wave propagation speed 𝑐 in positive and negative 𝑥−direction, respectively.

For a unique solution of the functions 𝑓 and 𝑔, the initial conditions have to be evaluated. Let

$$u\_0(\mathbf{x}) = u(\mathbf{x}, 0), \ v\_0(\mathbf{x}) = \left. \frac{\partial u(\mathbf{x}, t)}{\partial t} \right|\_{t=0} \tag{2.25}$$

be the initial displacement and velocity for −∞ ≤ 𝑥 ≤ ∞ at 𝑡 <sup>=</sup> 0. Then, the initial value problem is solved with [72]

$$u(\mathbf{x},t) = \frac{1}{2} \left[ u\_0(\mathbf{x} - ct) + u\_0(\mathbf{x} + ct) + \frac{1}{c} \int\_{\mathbf{x} - ct}^{\mathbf{x} + ct} v\_0(\mathbf{\tilde{x}}) \, d\mathbf{\tilde{x}} \right]. \tag{2.26}$$

The initial displacement 𝑢0(𝑥) splits up in two parts of half magnitude each of which are traveling in positive and negative 𝑥−direction. Moreover, due to the finite wave propagation speed, the initial velocity at 𝑥 <sup>=</sup> 𝑥¯ only influences the displacements in the interval 𝑥 ∈ [𝑥¯ <sup>−</sup> 𝑐𝑡, 𝑥¯ <sup>+</sup> 𝑐𝑡].

### **2.3.2 Transition and Boundary Conditions**

For practical applications, transition and boundary conditions have to be considered. An overview of the most common transition and boundary conditions is given in Fig. 2.3.

### **Transition Conditions**

The aim is to find a coherence between the known incident force wave 𝑁i(𝑥tr, 𝑡) and its reflected 𝑁r(𝑥tr, 𝑡) and transmitted parts 𝑁t(𝑥tr, 𝑡). Therefore, an incident wave

**Figure 2.3:** Two homogeneous elements with a fixed end at <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup>1, a free end at <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup><sup>2</sup> and the transition at 𝑥 <sup>=</sup> 𝑥tr.

𝑓i(𝑥 <sup>−</sup> 𝑐1𝑡), which travels in the first element toward the transition, is considered. At the transition, the incident wave is partly reflected and transmitted. This leads to traveling waves 𝑔r(𝑥 <sup>+</sup> 𝑐1𝑡) and 𝑓t(𝑥 <sup>−</sup> 𝑐2𝑡) in the first respectively second element. In general, the displacements of the two elements are described by

$$u\_k(\mathbf{x}, t) = f\_k(\mathbf{x} - c\_k t) + g\_k(\mathbf{x} + c\_k t), \ k = 1, 2,\tag{2.27}$$

which simplifies for the given consideration to

$$u\_1(\mathbf{x}, t) = f\_1(\mathbf{x} - c\_1 t) + g\_1(\mathbf{x} + c\_1 t),\tag{2.28}$$

$$u\_2(\mathbf{x}, t) = f\_t(\mathbf{x} - c\_2 t). \tag{2.29}$$

Moreover, two transition conditions have to be stated which are continuity of displacement

$$
\mu\_1\left(\mathbf{x}\_{\text{tr}}, t\right) = \mu\_2\left(\mathbf{x}\_{\text{tr}}, t\right) \tag{2.30}
$$

and force equilibrium

$$N\_{\rm l}(\mathbf{x}\_{\rm tr}, t) + N\_{\rm r}(\mathbf{x}\_{\rm tr}, t) = N\_{\rm l}(\mathbf{x}\_{\rm tr}, t). \tag{2.31}$$

The displacement 𝑢1(𝑥tr, 𝑡) is composed of an incident part 𝑢i(𝑥tr, 𝑡) and a reflected one 𝑢r(𝑥tr, 𝑡), whereas 𝑢2(𝑥tr, 𝑡) only consists of the transmitted part 𝑢t(𝑥tr, 𝑡). In terms of traveling waves, Eq. (2.30) is written as

$$\underbrace{f\_{\mathbf{l}}\left(\mathbf{x}\_{\text{tr}} - c\_{1}t\right)}\_{\mathbf{u}\_{\text{l}}(\mathbf{x}\_{\text{tr}},t)} + \underbrace{g\_{\mathbf{r}}\left(\mathbf{x}\_{\text{tr}} + c\_{1}t\right)}\_{\mathbf{u}\_{\text{l}}(\mathbf{x}\_{\text{tr}},t)} = \underbrace{f\_{\mathbf{l}}\left(\mathbf{x}\_{\text{tr}} - c\_{2}t\right)}\_{\mathbf{u}\_{\text{l}}(\mathbf{x}\_{\text{tr}},t)}.\tag{2.32}$$

As the displacement at <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup>tr is continuous, this is also true for the velocity

$$\left. \frac{\partial u\_1(\mathbf{x}, t)}{\partial t} \right|\_{\mathbf{x} = \mathbf{x}\_{\text{tr}}} = \left. \frac{\partial u\_2(\mathbf{x}, t)}{\partial t} \right|\_{\mathbf{x} = \mathbf{x}\_{\text{tr}}} \text{ \tag{2.33}$$

$$-c\_1 f\_\mathbf{i}'(\mathbf{x}\_{\text{tr}} - c\_1 t) + c\_1 g\_\mathbf{r}'(\mathbf{x}\_{\text{tr}} + c\_1 t) = -c\_2 f\_\mathbf{i}'(\mathbf{x} - c\_2 t),\tag{2.34}$$

$$-c\_1 \underbrace{E\_1 A\_1 f\_{\mathbf{i}}'(\mathbf{x}\_{\mathrm{tr}} - c\_1 t)}\_{N\_\mathrm{I}(\mathbf{x}\_{\mathrm{tr}}, t)} + c\_1 \underbrace{E\_1 A\_1 g\_{\mathbf{r}}'(\mathbf{x}\_{\mathrm{tr}} + c\_1 t)}\_{N\_\mathrm{I}(\mathbf{x}\_{\mathrm{tr}}, t)} = -c\_2 \underbrace{\frac{E\_1 A\_1}{E\_2 A\_2}}\_{N\_\mathrm{I}(\mathbf{x}\_{\mathrm{tr}}, t)} \underbrace{E\_2 A\_2 f\_{\mathbf{i}}'(\mathbf{x}\_{\mathrm{tr}} - c\_2 t)}\_{N\_\mathrm{I}(\mathbf{x}\_{\mathrm{tr}}, t)}\tag{2.35}$$

from which

$$-N\_{\rm l}(\mathbf{x}\_{\rm tr},t) + N\_{\rm r}(\mathbf{x}\_{\rm tr},t) = -\frac{Z\_1}{Z\_2}N\_{\rm l}(\mathbf{x}\_{\rm tr},t) \tag{2.36}$$

follows. Equation (2.36), together with Eq. (2.31) results in the final expressions

$$\mathbf{N\_{i}(x\_{\rm tr},t)} = \frac{2Z\_{2}}{Z\_{2} + Z\_{1}} \mathbf{N\_{i}(x\_{\rm tr},t)} = T\_{1,2} \mathbf{N\_{i}(x\_{\rm tr},t)},\tag{2.37}$$

$$\mathbf{N\_{r}(x\_{\text{tr}},t)} = \frac{Z\_{2} - Z\_{1}}{Z\_{2} + Z\_{1}} \mathbf{N\_{i}(x\_{\text{tr}},t)} = \mathbf{R\_{1,2}} \mathbf{N\_{i}(x\_{\text{tr}},t)}\tag{2.38}$$

for the transmitted and reflected force waves depending on the incident force. The transmission and reflection factors from element 𝑖 to element 𝑗 are defined by

$$T\_{l\_{\parallel}} = \frac{2Z\_{\parallel}}{Z\_{\parallel} + Z\_{\parallel}},\tag{2.39}$$

$$\mathcal{R}\_{l,j} = \frac{\mathbf{Z}\_{\parallel} - \mathbf{Z}\_{l}}{\mathbf{Z}\_{\parallel} + \mathbf{Z}\_{l}}.\tag{2.40}$$

#### **Boundary Conditions**

Along with the transition conditions, the boundary conditions are necessary to determine the solution for <sup>𝑢</sup>(𝑥, 𝑡). In Fig. 2.3, the boundary conditions are a fixed end at <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup><sup>1</sup> and a free end at 𝑥 <sup>=</sup> 𝑥2.

At the fixed end, the displacement is identical zero, i. e.

$$u(\mathbf{x}\_1, t) \equiv \mathbf{0},\tag{2.41}$$

whereas at the free end, the force is identical zero, i. e.

$$N(x\_2, t) \equiv 0.\tag{2.42}$$

The corresponding reflected forces at the fixed respectively free end are obtained by taking the limit <sup>𝑍</sup><sup>2</sup> → ∞ (fixed end) or setting <sup>𝑍</sup><sup>2</sup> <sup>=</sup> <sup>0</sup> (free end) in Eq. (2.38)

$$\mathbf{N\_{r}(x\_{1},t)} = \mathbf{N\_{l}(x\_{1},t)},\tag{2.43}$$

$$N\_{\mathbf{f}}(\mathbf{x}\_2, t) = -N\_{\mathbf{i}}(\mathbf{x}\_2, t). \tag{2.44}$$

# **3 Impacting Rods: Direct Problem**

The aim of this chapter is to verify the wave propagation method. Initially, the wave propagation method is presented. Its main equations for one time step are derived. Possible applications of the wave propagation method are demonstrated by two examples. Afterward, the wave propagation method is summarized. Subsequently, the modal expansion method is applied to a hit of two rods of constant and exponential increasing cross section area. After determining the eigenfunctions of the problem, the initial conditions are incorporated by using an orthogonality relation. At the end of the chapter, the results that are determined with both methods, are compared. It is shown that the wave propagation method provides more accurate results if a finite number of terms is used.

## **3.1 Wave Propagation Method**

The direct solution of the impacting rod problem has been discussed extensively in literature, e.g. in references [31, 38]. One approach is to model the rods with concentrated masses and stiffnesses. This also applies to the well known finite element method (FEM). The FEM delivers accurate results in a wide range of applications. However, the method reveals disadvantages when wave problems are considered: On the one hand, the accelerations are assumed to be finite so that they are numerically integrable. This is the reason why strong discontinuities in force and velocity, which often occur in impacting problems, are not displayed correctly. On the other hand, the FEM does not take into account that the wavefront is traveling through the elements with finite speed. Thus, displacements are calculated at axial positions which are not reached by the wavefront yet. These disadvantages vanish when the system is modeled by distributed masses and stiffnesses which are described by the wave equation. This approach has been studied by Shorr in several publications, e.g. in references [66, 67]. The following derivation of the wave propagation method is partly based on his publications.

The main idea in modeling is to approximate the impedances of the rods by piecewise constant step functions. In Fig. 3.1, the elements of the rods are depicted for the special case of constant material parameters <sup>𝜌</sup> and 𝐸. Therefore, the impedances only change due to the changing cross section area. This discretization is advantageous since analytical solutions exist on elements with constant impedance. Moreover, the wave propagation speed <sup>𝑐</sup>𝑗 is constant on each element <sup>𝑗</sup>. The element lengths <sup>ℓ</sup>𝑗 are adjusted so that the waves travel through every element during a constant time step <sup>Δ</sup>𝑡

$$
\hbar \ell\_f = \sqrt{\frac{E\_f}{\rho\_f}} \Delta t = c\_f \Delta t.\tag{3.1}
$$

Thus, all the traveling waves reach the boundaries of the corresponding elements at the same time.

If external forces or loads are applied, the following assumptions pertain:


**Figure 3.1:** Impact of two non-uniform rods and their modeling.

As long as the impact force is compressive, the rods are calculated as if they were stuck together. When the force tends to become tensile, the rods separate. Therefore, the rods are treated as two individual rods with free ends at their former hitting surfaces.

### **3.1.1 Equations for One Time Step**

During one time step, the waves are traveling from one element border to the next element border. At the end of each time step at 𝑡 <sup>=</sup> 𝑡0, all elements are in equilibrium (see Fig. 3.2) which means that both the normal force 𝑁 and the particle velocity 𝑣 is constant within the elements. The objective of this section is to derive the equations which are required to determine the new equilibrium state one time step later at <sup>𝑡</sup> <sup>=</sup> <sup>𝑡</sup><sup>0</sup> <sup>+</sup> <sup>Δ</sup>𝑡. To this end, the transition forces and velocities are determined first. This enables the calculation of the incident force and velocity waves starting from the element borders. The new forces and velocities are finally obtained by a superposition of the incident waves and their old values.

**Figure 3.2:** Two elements in equilibrium state at the end of each time step.

In general, the forces and velocities of two adjacent elements are not equal at 𝑡 <sup>=</sup> 𝑡0, that is <sup>𝑁</sup>𝑗−<sup>1</sup> <sup>≠</sup> <sup>𝑁</sup>𝑗 and <sup>𝑣</sup>𝑗−<sup>1</sup> <sup>≠</sup> <sup>𝑣</sup>𝑗 . These inequalities cause waves to propagate into the elements starting from the transition (see Fig. 3.3a). The superscripts indicate that 𝑁 or 𝑣 is either acting on the right end (superscript <sup>+</sup>) or left end (superscript <sup>−</sup>) of the corresponding element.

**Figure 3.3:** Free body diagrams at the transition.

At the element transitions, force equilibrium and continuity is claimed. As a consequence of continuity, the velocities of the related element particles have to be equal as well. Mathematically, the claims are expressed by

$$\rm N\_{/-1}^{+} - N\_{/}^{-} = F\_{/ \prime}^{\*} \tag{3.2}$$

$$\boldsymbol{v}\_{/-1}^{'} = \boldsymbol{v}\_{/}^{'} \tag{3.3}$$

with the external force 𝐹 <sup>∗</sup> and the transition forces 𝑁<sup>+</sup> 𝑗−1 , 𝑁− 𝑗 and velocities 𝑣 + 𝑗−1 , 𝑣− 𝑗 acting on the element borders. During the travel time ¯<sup>𝑡</sup> ∈ [𝑡0, 𝑡<sup>0</sup> <sup>+</sup> <sup>Δ</sup>𝑡], the waves have propagated the distances <sup>𝑐</sup>𝑗−<sup>1</sup> ¯<sup>𝑡</sup> and <sup>𝑐</sup>𝑗 ¯𝑡, respectively. In Fig. 3.3b, the free body diagrams of the parts that are covered by the incident waves are depicted. Applying the momentum conservation law on the free body diagram of element 𝑗

$$\begin{aligned} \left(\mathbf{N}\_{\rangle} - \mathbf{N}\_{\rangle}^{-}\right) \overline{t} &= A\_{\rangle} \rho\_{\rangle} c\_{\rangle} \overline{t} \left(v\_{\rangle}^{-} - v\_{\rangle}\right), \\ \mathbf{N}\_{\rangle}^{-} &= \mathbf{N}\_{\rangle} - Z\_{\rangle} \left(v\_{\rangle}^{-} - v\_{\rangle}\right) \end{aligned} \tag{3.4}$$

leads to an equation for the transition force acting on the left end of element 𝑗. Analogously, the transition force acting on the right end of element 𝑗 <sup>−</sup> <sup>1</sup> is derived

$$N\_{j-1}^{+} = N\_{j-1} + Z\_{j-1} \left( v\_{j-1}^{+} - v\_{j-1} \right) . \tag{3.5}$$

By evaluating the transition conditions (Eqs. (3.2), (3.3 )) together with the above derived equations (Eqs. (3.4), (3.5)), the unknown transition forces 𝑁<sup>+</sup> 𝑗 , 𝑁− 𝑗 and velocities 𝑣 + 𝑗 , 𝑣− 𝑗

$$N\_{j}^{+} = \frac{\left(\left(\upsilon\_{j+1} - \upsilon\_{j}\right)Z\_{j+1} + N\_{j+1} + F\_{j+1}^{\*}\right)Z\_{j} + N\_{j}Z\_{j+1}}{Z\_{j+1} + Z\_{j}},\tag{3.6}$$

$$Z\_{j} = \upsilon\_{j+1}Z\_{j+1} - \upsilon\_{j}Z\_{j+1} + \upsilon\_{j}Z\_{j+1} + F\_{j+1}^{\*}$$

$$v\_{\uparrow}^{+} = \frac{Z\_{\uparrow}v\_{\uparrow} + Z\_{\uparrow+1}v\_{\uparrow+1} - N\_{\uparrow} + N\_{\uparrow+1} + F\_{\uparrow+1}^{\*}}{Z\_{\uparrow+1} + Z\_{\uparrow}},\tag{3.7}$$

$$N\_f^- = \frac{\left(\left(\upsilon\_f - \upsilon\_{f-1}\right)Z\_{f-1} + N\_{f-1} - F\_f^\*\right)Z\_f + N\_f Z\_{f-1}}{Z\_{f-1} + Z\_f},\tag{3.8}$$

$$Z\_{f-1} + Z\_{f-1} = \upsilon\_{f-1}Z\_{f-1} + N\_f Z\_{f-1} + F\_f^\*$$

$$v\_{j}^{-} = \frac{Z\_{j}v\_{j} + Z\_{j-1}v\_{j-1} + N\_{j} - N\_{j-1} + F\_{j}^{\*}}{Z\_{j-1} + Z\_{j}} \tag{3.9}$$

are determined after the subscripts have been manipulated. The force difference <sup>Δ</sup>𝑁 between the transition forces 𝑁<sup>+</sup> 𝑗 , 𝑁− 𝑗 and the old force <sup>𝑁</sup>𝑗 leads to traveling force waves

$$
\Delta \mathbf{N}\_{\text{f,L}} = \mathbf{N}\_{\text{f}}^{+} - \mathbf{N}\_{\text{f}}.\tag{3.10}
$$

$$
\Delta \mathbf{N}\_{\rangle, \mathbf{R}} = \mathbf{N}\_{\rangle}^{-} - \mathbf{N}\_{\rangle} \tag{3.11}
$$

as can be seen in Fig. 3.4. The subscripts denote if the waves are traveling to the left

**Figure 3.4:** Incident waves in element 𝑗 and in its adjacent elements.

(subscript L) or right (subscript R). Analogously, the difference in velocity <sup>Δ</sup>𝑣 leads to traveling velocity waves

$$
\Delta v\_{\rangle,\mathcal{L}} = v\_{\rangle}^{+} - v\_{\rangle,\mathcal{L}} \tag{3.12}
$$

$$
\Delta \boldsymbol{\upsilon}\_{/\mathcal{R}} = \boldsymbol{\upsilon}\_{/}^{-} - \boldsymbol{\upsilon}\_{/}.\tag{3.13}
$$

The new force <sup>𝑁</sup>𝑗,new and velocity <sup>𝑣</sup>𝑗,new are finally obtained by a superposition of the incident waves and the old values

$$N\_{\rm/new} = N\_{\rm f} + \Delta N\_{\rm f,L} + \Delta N\_{\rm f,R} \tag{3.14}$$

$$
\sigma\_{\text{j,new}} = \upsilon\_{\text{j}} + \Delta \upsilon\_{\text{j,L}} + \Delta \upsilon\_{\text{j,R}}.\tag{3.15}
$$

### **3.1.2 Boundary Conditions**

The most common boundary conditions (given boundary force, given velocity/displacement, energy sink) are presented.

For **given boundary forces** 𝑁<sup>−</sup> 𝑛 , 𝑁+ 𝑛˜ , the subsequent equations for the velocities are derived by manipulating Eqs. (3.4) and (3.5)

$$
\upsilon\_n^- = \frac{N\_n - N\_n^-}{Z\_n} + \upsilon\_{n\prime} \tag{3.16}
$$

$$
\upsilon\_{\tilde{n}}^{+} = \frac{N\_{\tilde{n}}^{+} - N\_{\tilde{n}}}{Z\_{\tilde{n}}} + \upsilon\_{\tilde{n}}.\tag{3.17}
$$

For **given velocities** 𝑣 − 𝑛 , 𝑣+ 𝑛˜ , the resulting forces are determined by rearranging Eqs. (3.16) and (3.17)

$$\begin{aligned} \mathbf{N}\_n^- &= \mathbf{N}\_n + Z\_n \left(\boldsymbol{\upsilon}\_n - \boldsymbol{\upsilon}\_n^-\right), \\ \dots, \dots, \dots &= \dots \dots \dots \end{aligned} \tag{3.18}$$

$$N\_{\vec{\mathbb{R}}}^{+} = N\_{\vec{\mathbb{R}}} + Z\_{\vec{\mathbb{R}}} \left( \upsilon\_{\vec{\mathbb{R}}}^{+} - \upsilon\_{\vec{\mathbb{R}}} \right) . \tag{3.19}$$

For a **given energy sink**, i.e. that no reflecting waves occur at the boundaries, the equations are

$$N\_{\underline{n}}^{-} = N\_{\underline{n}\prime} \qquad N\_{\underline{n}}^{+} = N\_{\tilde{\underline{n}}}.\tag{3.20}$$

The corresponding velocities are obtained by substituting the determined forces into Eqs. (3.16) and (3.17).

### **3.1.3 Examples**

The range of applications of the wave propagation method is demonstrated by applying the method to two examples. The first example is a textbook example with one jump in impedance. The second example is more complex, including variations in geometry and material parameters. Moreover, friction is considered and its influence on the simulation results is analyzed.

#### **First Example**

The example presented in Fig. 3.5 is investigated. Rod 1 with initial velocity <sup>𝑣</sup>in is hitting Rod 2 which is initially at rest. The impedance function of Rod 2 is constant (𝑍˜ <sup>=</sup> 𝑍0). Rod 1 consists of two parts with impedance <sup>𝑍</sup><sup>0</sup> and <sup>𝑍</sup>1, respectively. The simulation is run with the parameters listed in Tab. 3.1. The corresponding results of the first

**Figure 3.5:** Setup of the first example. All dimensions in m.


**Table 3.1:** Parameters used for the first example.

example are presented in Fig. 3.6. The displacement (left) and the normal force (right) are depicted at three time instants. The transition between Rod 1 and Rod 2 is indicated by a vertical dashed line.

At the first time instant, compressive force waves are traveling in both rods. The part of the wavefront in Rod 1, which has been transmitted at the impedance jump of Rod 1, is traveling to the left. The reflected part is moving toward the transition between Rod 1 and Rod 2. The waves that have not been covered by the wavefronts yet are either traveling with initial velocity in Rod 1 or are at rest in Rod 2.

At the second time instant, the part of the incident wave in Rod 1, that has been reflected at the impedance jump, decreases the absolute value of the normal force in Rod 2 when it reaches the transition zone. Shortly after, when the incident wave in Rod 1 has traveled back and forth along Rod 1, the rods separate. Thus, both rods move independently which leads to discontinuous displacements of the rods at 𝑥 <sup>=</sup> 0. After the separation, the normal forces periodically change in both rods.

At the third time instant, the incident waves in Rod 2 have been reflected at its free end. Therefore, the compressive force waves traveling to the right have been converted to tensile force waves traveling to the left. The tensile force waves are pulling the particles in positive 𝑥−direction. The displacements in Rod 2 are lowest in the region that has not been reached by the reflected tensile force waves yet. Before the compressive waves have been reflected at the free end, the displacements in this region have been largest (see second time instant). This procedure is permanently repeated and leads to continuously increasing displacements in Rod 2. Physically, the results are comprehensible as Rod 2 is moving in positive 𝑥−direction after the impact. The displacements and normal forces in Rod 1 may be interpreted analogously.

**Figure 3.6:** Results of the first example. The displacements (left) and the forces (right) are presented at 𝑡 <sup>=</sup> <sup>0</sup>.064 ms (first row), 𝑡 <sup>=</sup> <sup>0</sup>.259 ms (second row) and 𝑡 <sup>=</sup> <sup>0</sup>.549 ms (third row).

#### **Complex Example**

The setup depicted in Fig. 3.7 is investigated where Rod 1 with initial velocity <sup>𝑣</sup>in hits Rod 2 with initial velocity <sup>𝑣</sup>˜in. Both rods are supported by two plain bearings (see Fig. 3.7a). It is assumed that Rod 1 can move freely, whereas friction is considered at the plain bearings of Rod 2. The first plain bearing of Rod 2 is located between <sup>𝑥</sup><sup>2</sup> and <sup>𝑥</sup>3, the second is located between <sup>𝑥</sup><sup>4</sup> and <sup>𝑥</sup>5. With the mass

$$m = \sum\_{\{\rho = 1\}}^{\tilde{n}} \tilde{A}\_{\{\tilde{\rho}\}} \tilde{\ell}\_{\{\rho\}} \tag{3.21}$$

of Rod 2, the total friction force

$$F\_{\mathbb{R}} = \mu m \,\overline{\text{g}}\tag{3.22}$$

is calculated, where <sup>𝜇</sup> is the friction coefficient and 𝑔¯ is the gravitational constant. If 𝑛<sup>ˆ</sup> is the total number of nodes 𝑘 in contact with the plain bearings of Rod 2, then it is assumed that the total friction force is uniformly distributed

$$F\_k^\* = -\frac{F\_\mathcal{R}}{\hat{\mu}} \operatorname{sgn} \left( \frac{1}{2} (\tilde{v}\_k + \tilde{v}\_{k-1}) \right). \tag{3.23}$$

The cross section area of Rod 1 is constant (𝐴 <sup>=</sup> 𝐴0) as can be seen in Fig. 3.7a. The cross section area of Rod 2 linearly increases from <sup>𝐴</sup><sup>0</sup> to <sup>𝐴</sup><sup>1</sup> between <sup>𝑥</sup><sup>1</sup> and <sup>𝑥</sup>2. From <sup>𝑥</sup><sup>2</sup> to <sup>𝑥</sup>3, the cross section area is constant (<sup>𝐴</sup> <sup>=</sup> <sup>𝐴</sup>1). Between <sup>𝑥</sup><sup>3</sup> and <sup>𝑥</sup>4, <sup>𝐴</sup> is caluclated as a function of 𝑥

$$A(\mathbf{x}) = A\_1 - 0.3A\_0 \sin \left( \pi \frac{\mathbf{x} - \mathbf{x}\_3}{\mathbf{x}\_4 - \mathbf{x}\_3} \right). \tag{3.24}$$

Subsequently, there is a constant section (<sup>𝐴</sup> <sup>=</sup> <sup>𝐴</sup>1) between <sup>𝑥</sup><sup>4</sup> and <sup>𝑥</sup><sup>5</sup> followed by a linearly decreasing part from <sup>𝐴</sup><sup>1</sup> to <sup>𝐴</sup><sup>0</sup> between <sup>𝑥</sup><sup>5</sup> and <sup>𝑥</sup>6. In addition to the geometry variations, the material paramaters vary between <sup>𝑥</sup><sup>3</sup> and <sup>𝑥</sup><sup>4</sup> (see. Fig. 3.7b). This section may be interpreted as a rod part of different material (𝐸 <sup>=</sup> <sup>1</sup>.<sup>5</sup> 𝐸0, <sup>𝜌</sup> <sup>=</sup> <sup>1</sup>.<sup>2</sup> <sup>𝜌</sup>0) that is screwed together with the other rod parts of Rod 2. For example in rotary drilling, heavy weight drill rods are applied to reduce fatigue [16]. The geometry variations together with the material parameter variations lead to the impedance function depicted in Fig. 3.7c. Due to the changing material parameters in 𝑥 ∈ [𝑥3, 𝑥4], the impedance function increases in this region which leads to jumps at <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup><sup>3</sup> and <sup>𝑥</sup> <sup>=</sup> <sup>𝑥</sup>4.

**Figure 3.7:** Setup of the complex example.


**Table 3.2:** Parameters used for the complex example.

The simulation is run with the parameters listed in Tab. 3.2. In Fig. 3.8, the corresponding results for the displacement (left) and the normal force (right) are depicted at three time instants. The transition between Rod 1 and Rod 2 is indicated by a vertical dashed line. At the first time instant, the waves have propagated both into Rod 1 and Rod 2. The parts, that have not been covered by the wave front yet, travel with their initial velocities. At the second time instant, the wave front in Rod 1 is traveling toward the transition after it has been reflected at its free end. In Rod 2, two jumps in normal force may be observed. The first jump at 𝑥 <sup>≈</sup> <sup>0</sup>.<sup>14</sup> <sup>m</sup> results from the reflection of the incident wave front at the impedance jump at 𝑥 <sup>=</sup> 𝑥3. The transmitted part of the incident wave front is traveling at <sup>𝑥</sup> <sup>≈</sup> <sup>0</sup>.<sup>79</sup> m. Since the wave propagation speed is larger between <sup>𝑥</sup><sup>3</sup> and <sup>𝑥</sup><sup>4</sup> compared to the other parts, the transmitted part of the incident wave front has traveled further than the reflected part. At the last time instant, Rod 1 has separated from Rod 2. After the separation, free end boundary conditions are applied at the former hitting surface at 𝑥 <sup>=</sup> 0. Therefore, the normal force at 𝑥 <sup>=</sup> <sup>0</sup> is zero for both rods. Moreover, the rods can move independently in negative 𝑥−direction which explains the displacement jump at 𝑥 <sup>=</sup> 0. The particles of Rod 2 at axial position 𝑥 <sup>&</sup>gt; <sup>1</sup>.<sup>38</sup> <sup>m</sup> have not been reached by the wave front which is why the absolute value of the displacement is maximum and the normal forces are close to zero in this region.

The normal forces are not exactly zero in this region, since the friction forces of the second bearing position have already influenced the normal forces. However, the influence of the friction forces on the results is very little. The reason is that the impact force is much larger compared to the friction forces. For the given example, the absolute value of the initial impact force is approximately 2250 times the absolute value of the total friction force. It can be concluded that if short travel times are considered, the friction forces are

**Figure 3.8:** Results of the complex example. The displacements (left) and the forces (right) are presented at 𝑡 <sup>=</sup> <sup>0</sup>.070 ms (first row), 𝑡 <sup>=</sup> <sup>0</sup>.146 ms (second row) and 𝑡 <sup>=</sup> <sup>0</sup>.255 ms (third row).

negligible. In the long run, however, the friction forces significantly decrease the total energy <sup>𝐸</sup>tot of the system as can be seen in Fig. 3.9. The total energy

$$\begin{split} E\_{\text{tot}} &= \frac{1}{2} \sum\_{j=1}^{n} A\_{j} \rho\_{j} \ell\_{j} \upsilon\_{j}^{2} + \frac{1}{2} \sum\_{j=1}^{n} \frac{N\_{j}^{2}}{E\_{j} A\_{j}} \ell\_{j} + \frac{1}{2} \sum\_{j=1}^{n} \tilde{A}\_{j} \tilde{\rho}\_{j} \tilde{\ell}\_{j} \tilde{\upsilon}\_{j}^{2} + \frac{1}{2} \sum\_{j=1}^{n} \frac{\tilde{N}\_{j}^{2}}{\tilde{E}\_{j} \tilde{A}\_{j}} \tilde{\ell}\_{j} \\ &= E\_{\text{kin,Rad}} \, 1 + E\_{\text{pot,Rad}} \, 1 + E\_{\text{kin,Rad}} \, 2 + E\_{\text{pot,Rad}} \, 2 \end{split} \tag{3.25}$$

consists of the kinetic and potential energy of both rods. For the evaluation, 100 iterations have been considered which means that the wave has traveled 100 times forth and back in Rod 2. Without friction, the system is conservative and the total energy remains constant. With friction, the total energy linearly decreases from its maximum value at the beginning to approximately 128 J at the end. Thus, in the long run, friction has to be considered.

The examples show that the presented wave propagation method is applicable to the longitudinal impact of rods with arbitrary material parameters and geometries. In addition, external forces that enable to model e.g. friction may be considered. The wave propagation method will be used as the basis for the recursive method, which is presented in the next chapter.

**Figure 3.9:** Comparison between the total energy <sup>𝐸</sup>tot with and without the consideration of friction in Rod 2. One iteration means that the wave has traveled back and forth in Rod 2.

### **3.1.4 Summary**

In Fig. 3.10, the summary of the method is depicted. The chart is divided into three parts (application, preprocessing, simulation).

$$\begin{aligned} t &= 0 \quad : \quad & N\_{f\prime}, v\_f \\ & & F\_{f\prime}^\* N\_{f\prime}^+, N\_{f\prime}^-, v\_f^+, v\_f^- + \text{BC} \\ t &= \Delta t \quad : \quad & N\_{f, \text{new}}, v\_{f, \text{new}} \end{aligned}$$

**Figure 3.10:** Summary of the wave propagation method.

A possible application of the presented method is the impact of non-uniform rods at initial velocities <sup>𝑣</sup>in and <sup>𝑣</sup>˜in, respectively. The rods are subdivided into piecewise constant impedances in the preprocessing step. Moreover, the element lengths are adjusted so that it takes one time step to travel through each element. Based on the initial displacements and velocities, the initial forces <sup>𝑁</sup>𝑗 and velocities <sup>𝑣</sup>𝑗 at <sup>𝑡</sup> <sup>=</sup> <sup>0</sup> are determined. Afterward, the simulation starts with setting the current external forces 𝐹 ∗ 𝑗 . Furthermore, the transition forces 𝑁<sup>+</sup> 𝑗 , 𝑁− 𝑗 and velocities 𝑣 + 𝑗 , 𝑣− 𝑗 are determined and the current boundary conditions are set. The time step ends by calculating the new forces <sup>𝑁</sup>𝑗,new and velocities <sup>𝑣</sup>𝑗,new in equilibrium. This procedure is repeated until the desired simulation time is achieved.

The computational effort of the method is 𝒪(𝑛 ∗ ) with 𝑛 <sup>∗</sup> being the total number of rod elements. This enables fast calculations using many elements.

## **3.2 Modal Expansion Method**

As discussed in Sec. 2.2, analytical solutions exist for the longitudinal wave propagation in rods with varying cross section area. However, the possible cross sectional area variations are very limited. The aim of this section is to apply the modal expansion method to the problem given in Fig. 3.11. Thus, a reference solution is generated to

**Figure 3.11:** Investigated setup.

verify the wave propagation method. The setup consists of two homogeneous rods. Rod 1 is hitting Rod 2 with initial velocity 𝑣in. The cross section area of Rod 2 increases exponentially and the displacements of the rods are denoted by 𝑢1(𝑥1, 𝑡) and 𝑢2(𝑥2, 𝑡), respectively.

#### **Eigenfunctions**

The eigenfunctions <sup>𝑈</sup>𝑛(𝑥𝑛) of the approach

$$\begin{split} u\_n(\mathbf{x}\_n, t) &= a\_n(\mathbf{x}) \text{ e}^{l(k\_n \mathbf{x}\_n - l\omega t)} \\ &= a\_n(\mathbf{x}) \text{ e}^{l\mathbf{k}\_{\text{fl}} \mathbf{x}\_n} \text{ e}^{-l\omega t} = \mathcal{U}\_n(\mathbf{x}\_n) \text{ e}^{-l\omega t}, \quad n = 1, 2 \end{split} \tag{3.26}$$

are indicated by substituting the related wave numbers and amplitudes of Tab. 2.1 into <sup>𝑈</sup>𝑛(𝑥𝑛) . This leads to the eigenfunctions

$$\mathcal{U}\_1(\mathbf{x}\_1) = \mathbb{C}\_1 \cos \left(\frac{\omega}{c\_0} \mathbf{x}\_1\right) + \mathbb{C}\_2 \sin \left(\frac{\omega}{c\_0} \mathbf{x}\_1\right),$$

$$\mathcal{U}\_2(\mathbf{x}\_2) = e^{\frac{-\xi\tau\_2}{2}} \left[ \mathbf{C}\_3 \cos \left( \sqrt{\left(\frac{\omega}{c\_0}\right)^2 - \frac{\xi^2}{4}} \mathbf{x}\_2 \right) + \mathbf{C}\_4 \sin \left( \sqrt{\left(\frac{\omega}{c\_0}\right)^2 - \frac{\xi^2}{4}} \mathbf{x}\_2 \right) \right] \tag{3.28}$$

of Rod 1 and Rod 2. The unknown coefficients <sup>𝐶</sup><sup>1</sup> . . . 𝐶<sup>4</sup> are determined by evaluating the boundary and transition conditions. At the transition, displacement and force continuity are claimed. Moreover, the left end of Rod 1 is free whereas the right end of Rod 2 is fixed. The resulting equations

$$\left.\frac{\mathrm{d}LI\_1}{\mathrm{d}x\_1}\right|\_{x\_1=0} = 0,\tag{3.29}$$

$$\mathcal{U}\_1(\ell\_1) = \mathcal{U}\_2(0),$$
 
$$\begin{array}{c} \dots \ \ \ \ \ \end{array} \tag{3.30}$$

$$\left. \frac{\mathbf{d}LU\_1}{\mathbf{d}x\_1} \right|\_{\mathbf{x}\_1=\ell\_1} = \left. \frac{\mathbf{d}U\_2}{\mathbf{d}x\_2} \right|\_{\mathbf{x}\_2=0} \, \tag{3.31}$$

$$\mathcal{U}\_2(\ell\_2) = 0\tag{3.32}$$

are written as a system of linear equations

$$
\begin{pmatrix} 0 & \alpha & 0 & 0 \\\\ \cos(\alpha \ell\_1) & \sin(\alpha \ell\_1) & -1 & 0 \\\\ -\alpha \sin(\alpha \ell\_1) & \alpha \cos(\alpha \ell\_1) & \frac{\zeta}{2} & -\beta \\\\ 0 & 0 & \cos(\beta \ell\_2) & \sin(\beta \ell\_2) \end{pmatrix} \begin{pmatrix} C\_1 \\ C\_2 \\ C\_3 \\ C\_4 \\ C\_4 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \tag{3.33}
$$

with

$$
\alpha = \frac{\alpha}{c\_0}, \quad \beta = \sqrt{\alpha^2 - \left(\frac{\xi}{2}\right)^2}. \tag{3.34}
$$

Evaluating the characteristic equation

$$a\left(2\cos(a\ell\_1)\cos(\beta\ell\_2)\beta + \cos(a\ell\_1)\sin(\beta\ell\_2)\xi - 2a\sin(a\ell\_1)\sin(\beta\ell\_2)\right) = 0\tag{3.35}$$

determines the infinite number of eigenfrequencies <sup>𝜔</sup>𝑘 from which their corresponding coefficients

$$\mathbf{C}\_{1,k} = \mathbf{C}\_{3,k} \frac{1}{\cos\left(\frac{\omega\_k}{c\_0}\right)},\tag{3.36}$$

$$\mathbf{C}\_{2,k} = \mathbf{0}\_{\prime} \tag{3.37}$$

$$\mathbf{C}\_{4,k} = -\mathbf{C}\_{3,k} \frac{1}{\tan\left(\left(\frac{\omega\_0}{c\_0}\right)^2 - \frac{1}{16}\right)}\tag{3.38}$$

are obtained. Thus, the eigenfunctions are determined and the approach is rewritten

$$\begin{split} \left\{ \begin{aligned} \mu\_{1}(\mathbf{x}\_{1},t) \\ \mu\_{2}(\mathbf{x}\_{2},t) \end{aligned} \right\} &= \sum\_{k=1}^{\infty} \begin{Bmatrix} \mathcal{U}\_{1,k}(\mathbf{x}\_{1}) \\ \mathcal{U}\_{2,k}(\mathbf{x}\_{2}) \end{Bmatrix} [\mathcal{C}\_{3,k}\tilde{\mathcal{C}}\_{5,k}\sin(\omega\_{k}t) + \mathcal{C}\_{3,k}\tilde{\mathcal{C}}\_{6,k}\cos(\omega\_{k}t)] \\ &= \sum\_{k=1}^{\infty} \begin{Bmatrix} \mathcal{U}\_{1,k}(\mathbf{x}\_{1}) \\ \mathcal{U}\_{2,k}(\mathbf{x}\_{2}) \end{Bmatrix} [\mathcal{C}\_{5,k}\sin(\omega\_{k}t) + \mathcal{C}\_{6,k}\cos(\omega\_{k}t)]. \end{split} \tag{3.39}$$

The remaining coefficients <sup>𝐶</sup>5,𝑘 , 𝐶6,𝑘 are obtained by evaluating the initial conditions.

#### **Initial Conditions and Orthogonality Relation**

For a PDE of the form

$$\mathcal{S}(\mathbf{x})\frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} + \mathcal{K}[u(\mathbf{x},t)] = 0, \ \mathbf{x} \in [0, \ell] \tag{3.40}$$

with linear differential operator 𝒦[·], Hagedorn and DasGupta [36] have derived the orthogonality relation

$$\int\_{0}^{l} \mathfrak{G}(\mathbf{x}) \mathcal{U}\_{l}(\mathbf{x}) \mathcal{U}\_{k}(\mathbf{x}) \, \mathrm{d}x = B\_{k} \delta\_{/k} \, \mathrm{d}x \tag{3.41}$$

where <sup>𝐵</sup>𝑘 is a constant. Applying this to the given problem of two impacting rods with varying cross section area and constant material parameters leads to

$$\mathcal{H}\_{l}(\mathbf{x}\_{l}) = \rho A\_{l}(\mathbf{x}\_{l}), \quad \mathcal{H}\_{l}[\cdot] = -\frac{\partial}{\partial \mathbf{x}\_{l}} \left( EA\_{l}(\mathbf{x}\_{l}) \frac{\partial}{\partial \mathbf{x}\_{l}} \right), \quad i = 1, 2 \tag{3.42}$$

and the orthogonality relation is

$$\begin{aligned} \int\_{0}^{\ell\_1} A\_1(\mathbf{x}\_1) \mathcal{U}\_{1,/}(\mathbf{x}\_1) \mathcal{U}\_{1,k}(\mathbf{x}\_1) \, \mathrm{d}x\_1 \\ \int\_{0}^{\ell\_2} A\_2(\mathbf{x}\_2) \mathcal{U}\_{2,/}(\mathbf{x}\_2) \mathcal{U}\_{2,k}(\mathbf{x}\_2) \, \mathrm{d}x\_2 = \begin{cases} B\_{k\prime} & j=k \\ 0, & \text{else}. \end{cases} \end{aligned} \tag{3.43}$$

Finally, the initial conditions for the displacements

$$
\begin{Bmatrix} \mu\_1(\mathbf{x}\_1, \mathbf{0})\\ \mu\_2(\mathbf{x}\_2, \mathbf{0}) \end{Bmatrix} = \sum\_{k=1}^{\infty} \begin{Bmatrix} \mathcal{U}\_{1,k}(\mathbf{x}\_1) \\ \mathcal{U}\_{2,k}(\mathbf{x}\_2) \end{Bmatrix} \mathbf{C}\_{6,k} \tag{3.44}
$$

and for the velocities

$$\begin{Bmatrix} v\_1(\mathbf{x}\_1, \mathbf{0})\\ v\_2(\mathbf{x}\_2, \mathbf{0}) \end{Bmatrix} = \sum\_{k=1}^{\infty} \begin{Bmatrix} \mathcal{U}\_{1,k}(\mathbf{x}\_1) \\ \mathcal{U}\_{2,k}(\mathbf{x}\_2) \end{Bmatrix} \mathbf{C}\_{5,k} \omega\_k \tag{3.45}$$

are evaluated and the orthogonality relation is applied. Thus, the remaining coefficients

$$\mathbf{C}\_{5,k} = \frac{1}{B\_k \omega\_k} \left[ \int\_0^{\ell\_1} v\_1(\mathbf{x}\_1, 0) A\_1(\mathbf{x}\_1) U\_{1,k}(\mathbf{x}\_1) \, d\mathbf{x}\_1 + \int\_0^{\ell\_2} v\_2(\mathbf{x}\_2, 0) A\_2(\mathbf{x}\_2) U\_{2,k}(\mathbf{x}\_2) \, d\mathbf{x}\_2 \right] \tag{3.46}$$

and

$$\mathbf{C}\_{6,k} = \frac{1}{B\_k} \left[ \int\_0^{l\_1} u\_1(\mathbf{x}\_1, 0) A\_1(\mathbf{x}\_1) \mathbf{l} I\_{1,k}(\mathbf{x}\_1) \, d\mathbf{x}\_1 + \int\_0^{l\_2} u\_2(\mathbf{x}\_2, 0) A\_2(\mathbf{x}\_2) \mathbf{l} I\_{2,k}(\mathbf{x}\_2) \, d\mathbf{x}\_2 \right] \tag{3.47}$$

are obtained.

## **3.3 Comparison between Wave Propagation and Modal Expansion Method**

For piecewise constant impedance functions of the rods, the wave propagation method yields exact solutions up to a certain time as a sum of a finite number of terms (closedform). As compared to this, the modal expansion method leads to approximate solutions for all time if a finite number of terms is used. Both methods are exact for all time if one makes use of an infinite number of terms. Therefore, if it is only of interest to determine the solution up to a certain time, the wave propagation method is more accurate. These statements are demonstrated by applying both methods to the example depicted in Fig. 3.11 and comparing their results. Rod 1 is hitting Rod 2 at the initial velocity 𝑣in. Moreover, the initial displacements 𝑢1(𝑥1, <sup>0</sup>) and 𝑢2(𝑥2, <sup>0</sup>) are zero. The other parameter values are listed in Tab. 3.3.


**Table 3.3:** Parameters of the example presented in Fig. 3.11.

Figure 3.12 shows the corresponding results of the displacements and forces at three different times. The hit takes place at 𝑥 <sup>=</sup> <sup>0</sup> which is highlighted by a vertical dashed line. To the left of the dashed line, the results of Rod 1 are plotted whereas the results of Rod 2 are depicted on the right side. The rods stick together as long as the force at the transition is compressional. This condition always applies to the presented time steps.

At the first time, the incident waves have propagated <sup>0</sup>.<sup>75</sup> <sup>m</sup> in both rods. The part of the rods that is further away than <sup>0</sup>.<sup>75</sup> <sup>m</sup> from the hitting surface is not influenced by the impact yet. Thus, these elements are either at rest (Rod 2) or traveling at initial velocity (Rod 1) as can be seen on the left graph of the first row in Fig. 3.12. Furthermore, there is no significant deviation between the displacement results. These curves match significantly better than the force curves. The reason is that the force curves are not continuous. In the zoom box, the overshooting of the modal expansion method result at the jump is highlighted. This behavior is characteristic for eigenfunction series and is well known in literature as the Gibbs phenomenon. The forces in the parts of the rods, that are not covered by their wavefronts yet, are still zero. However, the oscillations of the modal expansion method results suggest that the wave already traveled further than it actually did. Apart from the region close to the discontinuities, the force curves are in very good accordance.

The continuous increase in cross section area of Rod 2 leads to continuously increasing forces according to amount in Rod 2 and at the transition. Since the cross section area of Rod 1 is constant, the wave shape does not change while propagating along Rod 1. Thus, the impact force at 𝑡 <sup>=</sup> <sup>0</sup> conforms with the wavefront of Rod 1.

**Figure 3.12:** Comparison of results obtained with the modal expansion method versus results obtained with the wave propagation method. The displacements (left) and the forces (right) are presented at 𝑡 <sup>=</sup> <sup>0</sup>.<sup>145</sup> ms (first row), 𝑡 <sup>=</sup> <sup>0</sup>.290 ms (second row) and 𝑡 <sup>=</sup> <sup>0</sup>.435 ms (third row).

At the second time, the waves have been reflected at their respective boundaries and travel back toward the transition. The free end on the left side of Rod 1 leads to zero force at 𝑥 <sup>=</sup> <sup>−</sup>1, whereas the fixed end of Rod 2 results in zero displacement at 𝑥 <sup>=</sup> 1. When the wavefront of Rod 1 reaches its left end, the maximum displacement is achieved. Afterward, the reflected tensile force wave reduces the maximum deflection.

At the third time, the wavefronts have met at the transition and are traveling back to their ends. This time, the reflected waves yield lower force magnitudes in Rod 2 and greater magnitudes in Rod 1. The force jumps always occur at the curve bends of the displacements.

In summary, the comparison of the results obtained with the two methods shows very good accordance. The deviations at the region close to the force jumps are explainable with the Gibbs phenomenon. It is disadvantageous for the modal expansion method that a large number of eigenfunctions (2500 for the presented example) have to be considered to plot the jumps. Moreover, the modal expansion method will always yield approximate solutions as an infinite number of eigenfunctions can not be considered. It can be concluded that for impact problems which include force jumps, the wave propagation method is preferable. However, most importantly, the wave propagation method is verified by the modal expansion method.

# **4 Impacting Rods: Inverse Problem**

Inverse problems occur in many practical applications. A well-known example is the computer tomography where by detecting an attenuated X-Ray, the spatial X-Ray absorption is concluded. This enables to get non-destructive information about the scanned structure. The purpose of the presented method is to determine the impedance function of Rod 1 which leads to the desired force over time at an arbitrary position of Rod 2. To solve the inverse problem, the impedance function of Rod 2 and the initial and boundary conditions have to be known (see Fig. 4.1). For the solution of the direct problem, only the impedance functions of the rods and the initial and boundary conditions have to be given to calculate the force and displacement distribution in both rods.

**Figure 4.1:** Overview of given and calculated values of the direct and inverse problem.

The following chapter is structured as follows: Initially, the recursive method (RM) is presented. The term reursive refers to the way how the inverse problem is solved. The implementation of the main idea is presented and a condition of existence is derived. It is shown how the impedance changes correlate with the probability that a solution of the inverse problem exists. Moreover, it is shown that a prescribed force at an arbitrary position of the second rod may be transformed to the corresponding impact force. At the end of the chapter, several examples are presented, including the application in percussion drilling.

## **4.1 Recursive Method**

The RM is based on the wave propagation method presented in Sec. 3.1. Without loss of generality, the hit takes place at 𝑡 <sup>=</sup> <sup>0</sup> and Rod 1 is the rod with unknown impedance. In addition to the assumptions that are valid for the wave propagation method, some prerequisites have to be met to apply the method.

These **prerequisites** are:


For Rod 2, there are no prerequisites.

In Fig. 4.2a, a schematic sketch of the elements close to the impact zone is depicted. The rods hit each other at initial velocities <sup>𝑣</sup>in of Rod 1 and <sup>𝑣</sup>˜in of Rod 2. Moreover, the impedances of Rod 2 are given. The impedances of Rod 1 are determined according to the predefined impact force (target force) which is illustrated in Fig. 4.2b. The main idea of the method is to utilize that each element of Rod 1 influences the impact force <sup>𝑁</sup><sup>I</sup> at different moments for the first time. Immediately after the impact, the influence of the first element of Rod 1 can be observed in Fig. 4.2b. Two time steps later, the influence of the second element becomes obvious. If the first 𝑞 <sup>−</sup> <sup>1</sup> elements in Fig. 4.2 are set to perfectly match the target values, the question remains how the 𝑞th element has to be set to continuously match the target force. Initially, the simulation is carried out with the same impedance as the preceding element (𝑍𝑞 <sup>=</sup> <sup>𝑍</sup>𝑞−1). The difference between the target force <sup>𝑁</sup>𝑞 (dashed line) and the force <sup>𝑁</sup>¯ 𝑞 obtained with equal impedances is characterized by <sup>Δ</sup>𝑁𝑞 . The occurring force difference is balanced by an adjustment of the impedance <sup>𝑍</sup>𝑞 . The only possibility that the impedance <sup>𝑍</sup>𝑞 influences the impact force <sup>𝑁</sup><sup>I</sup> at the considered time is if the initial impact force <sup>𝑁</sup><sup>1</sup> is directly transmitted to the transition between element 𝑞 <sup>−</sup> <sup>1</sup> and 𝑞 where it is reflected and directly transmitted back until it reaches the hitting surface. The mathematical implementation of the idea is presented in the following section.

**Figure 4.2:** Idea of the recursive method.

As the waves have to travel back and forth at each element before they reach the hitting surface, it is only possible to change <sup>𝑁</sup><sup>I</sup> every second time step by an impedance change of Rod 1. Therefore, the target force, the boundary force and the initial conditions are assumed to change only every second time step. However, this is not a serious limitation as by doubling the number of elements, the same discretization accuracy is obtained.

### **4.1.1 Prescribed Impact Force**

Initially, the formula to calculate the impedance of the first element is derived. Subsequently, the formula to recursively determine all the other elements of Rod 1 is presented. There is a different treatment of the first and all the other elements of Rod 1. The impedance of the first element is set to determine the appropriate initial impact force whereas for all the other elements, the initial impact force is given.

#### **First element**

The initial impact force

$$N\_1 = \frac{2Z\_1}{Z\_1 + \tilde{Z}\_1} \frac{(\tilde{v}\_{\rm in} - v\_{\rm in})\tilde{Z}\_1}{2} = T\_{\tilde{1},1}N\_0 \tag{4.1}$$

is obtained by inserting the initial conditions for the first element of Rod 1 and Rod 2 in Eq. (3.6). It is composed of the transmission factor

$$T\_{\tilde{1},1} = \frac{2Z\_1}{Z\_1 + \tilde{Z}\_1} \tag{4.2}$$

from the first element of Rod 2 to the first element of Rod 1 (see Eq. (2.39)) times the impact force

$$N\_0 = -\frac{(v\_{\rm in} - \tilde{v}\_{\rm in})\tilde{Z}\_1}{2}.\tag{4.3}$$

Rearranging Eq. (4.1) and solving for the unknown impedance <sup>𝑍</sup><sup>1</sup> yields

$$Z\_1 = -\frac{N\_1 \tilde{Z}\_1}{N\_1 + (\upsilon\_{\rm in} - \tilde{\upsilon}\_{\rm in})\tilde{Z}\_1}.\tag{4.4}$$

#### **Other elements**

The remaining elements (second until 𝑛th element) are determined by implementing the idea that every element influences the impact force at different moments for the first time. In Fig. 4.2, the first 𝑞 <sup>−</sup> <sup>1</sup> impedances are set to perfectly match the target force. The impedance of the next element <sup>𝑞</sup> is determined by first setting <sup>𝑍</sup>𝑞 <sup>=</sup> <sup>𝑍</sup>𝑞−<sup>1</sup> and determining the impact force <sup>𝑁</sup>¯ <sup>𝑞</sup> . Subsequently, the force difference <sup>Δ</sup>𝑁𝑞 between target force <sup>𝑁</sup>𝑞 and <sup>𝑁</sup>¯ 𝑞 is balanced by

$$\begin{split} \Delta N\_{q} = N\_{q} - \tilde{N}\_{q} &= N\_{1} \, T\_{1,2} \, T\_{2,3} \dots \, T\_{q-2,q-1} \, R\_{q-1,q} \, T\_{q-1,q-2} \dots \, T\_{3,2} \, T\_{2,1} \, T\_{1,\tilde{1}} \\ &= N\_{0} \, T\_{\tilde{1},1} \, T\_{1,2} \, T\_{2,3} \dots \, T\_{q-2,q-1} \, T\_{q-1,q-2} \dots \, T\_{3,2} \, T\_{2,1} \, T\_{1,\tilde{1}} \, R\_{q-1,q} \\ &= \left( T\_{\tilde{1},1} \, T\_{1,\tilde{1}} \right) \left( T\_{1,2} \, T\_{2,1} \right) \dots \left( T\_{q-2,q-1} \, T\_{q-1,q-2} \right) \, R\_{q-1,q} \, N\_{0,q} \end{split} \tag{4.5}$$

where <sup>𝑇</sup>𝑖,𝑗 , 𝑇𝑗,𝑖 are the transmission factors from element <sup>𝑖</sup> to element <sup>𝑗</sup> respectively from element <sup>𝑗</sup> to element <sup>𝑖</sup> and <sup>𝑅</sup>𝑞−1,𝑞 is the reflection factor from element <sup>𝑞</sup> <sup>−</sup> <sup>1</sup> to element 𝑞. There is only one wave train that influences the impact force due to the impedance change of <sup>𝑍</sup>𝑞 : The initial impact force <sup>𝑁</sup><sup>1</sup> which is directly transmitted to the (<sup>𝑞</sup> <sup>−</sup> <sup>1</sup>)th element (𝑇1,<sup>2</sup> . . . 𝑇𝑞−2,𝑞−<sup>1</sup> ) where it is reflected at the transition to the <sup>𝑞</sup>th element (𝑅𝑞−1,𝑞 ) and transmitted back to the impact surface (𝑇𝑞−1,𝑞−<sup>2</sup> . . . 𝑇1,1˜). Since

$$\mathbf{T}\_{l,\parallel}\mathbf{T}\_{\parallel,l} = \frac{4\mathbf{Z}\_{l}\mathbf{Z}\_{\parallel}}{(\mathbf{Z}\_{l}+\mathbf{Z}\_{\parallel})^{2}} = \mathbf{K}\_{l,\parallel} = \mathbf{K}\_{\parallel,l} \tag{4.6}$$

holds, the force difference can be written in a compact form

$$\begin{split} \Delta N\_{q} &= K\_{\overline{1},1} \, K\_{1,2} \dots \, K\_{q-2,q-1} \, R\_{q-1,q} \mathbf{N}\_{0} \\ &= K\_{q} \mathbf{R}\_{q-1,q} \mathbf{N}\_{0} \\ &= K\_{q} \frac{\mathbf{Z}\_{q} - \mathbf{Z}\_{q-1}}{\mathbf{Z}\_{q} + \mathbf{Z}\_{q-1}} \mathbf{N}\_{0} \end{split} \tag{4.7}$$

where

$$K\_q = K\_{\tilde{1},1} \prod\_{r=1}^{q-2} K\_{r,r+1} = K\_{\tilde{1},1} \prod\_{r=1}^{q-2} \frac{4Z\_r Z\_{r+1}}{(Z\_r + Z\_{r+1})^2}. \tag{4.8}$$

Rearranging Eq. (4.7) and solving for <sup>𝑍</sup>𝑞 leads to the final equation for the unknown impedance

$$Z\_q = Z\_{q-1} \frac{K\_q}{K\_q} \frac{N\_0 + \Delta N\_q}{N\_0 - \Delta N\_q} \,. \tag{4.9}$$

After determining <sup>𝑍</sup>𝑞 , the procedure repeats for the next element until all elements of Rod 1 are calculated recursively.

### **4.1.2 Transformation of the Prescribed Force**

If the force <sup>𝑁</sup>𝑝 is prescribed at an arbitrary position <sup>𝑝</sup> of Rod 2, it is possible to transform <sup>𝑁</sup>𝑝 to the corresponding impact force <sup>𝑁</sup><sup>I</sup> (see Fig. 4.3a). The presented procedure makes use of the linearity of the problem. Thus, the impact force is reconstructed by a sum of weighted time-shifted responses to a rectangular excitation. The prerequisites of the RM and the piecewise constant modeling of the rod's impedances lead to forces that only vary every second time step. The prescribed force <sup>𝑁</sup>𝑝 (𝑡) at position <sup>𝑝</sup> is defined by its factors <sup>𝑌</sup>𝑗

$$N\_p\left(t\right) = Y\_{\left|f\right.}\ t\_{\rm dt} + 2(j-1)\Delta t \le t < t\_{\rm dt} + 2j\Delta t,\ j = 1\ldots n \tag{4.10}$$

as can be seen in Fig. 4.3b.

#### 4 Impacting Rods: Inverse Problem

**Figure 4.3:** Recursive determination of the factors <sup>𝑋</sup><sup>𝑗</sup> of the impact force <sup>𝑁</sup><sup>I</sup> .

The aim is to determine the factors <sup>𝑋</sup>𝑗

$$N\_{\mathbf{I}}(t) = X\_{\mathbf{j}}, \ \mathbf{2}(j-1)\boldsymbol{\Delta t} \le t < \mathbf{2}j\boldsymbol{\Delta t}, \ j = 1\ldots n \tag{4.11}$$

of the unknown force <sup>𝑁</sup><sup>I</sup> which generate <sup>𝑁</sup>𝑝. To this end, a rectangular function

$$N\_{\rm I}(t) = \mathcal{H}(t) - \mathcal{H}(t - 2\Delta t),\tag{4.12}$$

with ℋ (·) being the Heaviside function, is applied. Its response at position 𝑝

$$\mathcal{G}\left(t\right) = \mathcal{G}\_{\left(\prime\right)}\;t\_{\text{dt}} + \mathcal{Q}(j-1)\Delta t \le t < t\_{\text{dt}} + \mathcal{Q}j\Delta t,\; j = 1\ldots n \tag{4.13}$$

is called rectangular response. In Fig. 4.3c, the rectangular response together with the desired force is depicted. The rectangular response is weighted so that it matches the first step of <sup>𝑁</sup>𝑝 (see Fig. 4.3d). Therefore, the first factor is determined as

$$X\_1 = \frac{Y\_1}{G\_1}.\tag{4.14}$$

All the other time-shifted rectangular responses arrive later and are not able to influence the first step. In addition to the first rectangular function, a second rectangular function

$$N\_{\rm I}(t) = \mathcal{H}(t - 2\Delta t) - \mathcal{H}(t - 4\Delta t),\tag{4.15}$$

which is time-shifted by <sup>2</sup>Δ𝑡, is applied (see Fig. 4.3e). The second rectangular response is weighted to match the second step. This enables the calculation of the second factor

$$X\_2 = \frac{1}{G\_1} \left( Y\_2 - X\_1 G\_2 \right). \tag{4.16}$$

The procedure is repeated until all factors <sup>𝑋</sup>𝑗 are determined recursively. In general, the formula for the determination of the 𝑗th factor is

$$X\_{\!/\!} = \frac{1}{G\_1} \left( Y\_{\!/\!} - \sum\_{k=1}^{\langle \!/ -1 \!\!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/ \!/} \right) . \tag{4.17}$$

### **4.1.3 Condition of Existence**

It is not always possible to generate a solution of the inverse problem for any given force over time. Physical solutions only exist if the obtained impedances in Eq. (4.9) are positive, i. e.

$$Z\_{\theta} > 0.\tag{4.18}$$

Therefore, the following condition

$$\frac{K\_q \, N\_0 + \Delta N\_q}{K\_q \, N\_0 - \Delta N\_q} > 0\tag{4.19}$$

holds which is equivalent to

$$\left|\Delta \mathbf{N}\_q\right| < \left|\mathbf{K}\_q \cdot \mathbf{N}\_0\right|.\tag{4.20}$$

The force <sup>𝑁</sup><sup>0</sup> is problem specific but independent of the determined element. By rewriting the factor <sup>𝐾</sup>𝑞 as

$$K\_q = K\_{\tilde{1},1} \prod\_{r=1}^{q-2} K\_{r,r+1} = K\_{\tilde{1},1} \prod\_{r=1}^{q-2} \left[ 1 - \left( \frac{Z\_r - Z\_{r+1}}{Z\_r + Z\_{r+1}} \right)^2 \right] \le 1,\tag{4.21}$$

it becomes obvious that the absolute value of the force reserve 𝐾𝑞 𝑁<sup>0</sup>  decreases if the consecutive impedances change. This means that the larger changes in impedance are, the lower are the force reserves. Therefore, lower force differences <sup>Δ</sup>𝑁𝑞  may be compensated as condition (4.20) has to hold. Physically, the factor <sup>𝐾</sup>𝑞 may be interpreted as the fraction of transmitted energy through the impactor. The more impedance steps are passed and the higher the jumps are, the more internal reflections occur whose energy cannot be used to compensate the force differences.

In Fig. 4.4, the factors <sup>𝐾</sup>𝑟,𝑟+<sup>1</sup> and <sup>𝐾</sup>𝑛 are analyzed for a constant change in impedance, that is <sup>𝑍</sup>𝑟+<sup>1</sup> <sup>=</sup> 𝑚𝑍𝑟 . The factor

$$K\_{\rm n} = \prod\_{r=1}^{n} \left[ 1 - \left( \frac{Z\_r - Z\_{r+1}}{Z\_r + Z\_{r+1}} \right)^2 \right] = \left[ 1 - \left( \frac{1 - m}{1 + m} \right)^2 \right]^n \tag{4.22}$$

describes the fraction of total transmitted energy for 𝑛 transitions and

$$K\_{r,r+1} = 1 - \left(\frac{1-m}{1+m}\right)^2\tag{4.23}$$

describes the fraction of total transmitted energy at one transition. As the fraction of total transmitted energy at one transition is independent of whether the wave is first transmitted from the element with lower impedance to the element with larger impedance or reverse, <sup>𝐾</sup>𝑟,𝑟+1(𝑚) <sup>=</sup> <sup>𝐾</sup>𝑟,𝑟+<sup>1</sup> 1 𝑚 applies. Therefore, only values 𝑚 <sup>≥</sup> <sup>1</sup> are plotted. For constant impedances (𝑚 <sup>=</sup> 1), all the energy is transmitted (see Fig. 4.4a). The larger the changes in impedance are, the lower is the fraction of transmitted energy. In Fig. 4.4b, it is shown that large jumps in impedance lead to disproportionally large decreases in <sup>𝐾</sup>𝑛. For example, one jump in impedance with <sup>𝑚</sup> <sup>=</sup> <sup>2</sup> leads to larger

**Figure 4.4:** Analysis of the continuously decreasing force reserve.

decreases in <sup>𝐾</sup>𝑛 than <sup>100</sup> jumps with <sup>𝑚</sup> <sup>=</sup> <sup>1</sup>.<sup>05</sup> although the initial impedance value only doubles in the first case and approximately 132-folds in the second case.

### **4.1.4 Summary**

The presented method applies to many inverse impact problems. The position, where the force over time is given, can be chosen arbitrarily if the prerequisites are met. It is only possible to prescribe the force over time for a time span which is twice the travel time of the waves through Rod 1. Afterward, there are no elements left to shape the waves. By applying the RM, the two main questions concerning inverse problems are answered: Firstly, there is a solution of the inverse problem (**existence**) if the condition in Eq. (4.20) holds for every element of Rod 1. Moreover, in case the condition is violated, it is possible to assess the maximum length of Rod 1 for which a solution of the inverse problem still exists. Secondly, the solution is unique (**uniqueness**) in terms of impedance. Another important aspect that is frequently discussed in inverse problems is the stability of the results with respect to disturbed measurements. However, as the prescribed force in the presented impacting rod problem is desired and not measured, there are no deviations which could lead to instabilities.

The calculation time of the algorithm (see Tab. 4.1) is 𝒪(𝑛 ∗2 ) with 𝑛 <sup>∗</sup> being the total number of rod elements. On a standard laptop, the solution is obtained within 1 s if less than 400 elements are chosen.


**Table 4.1:** Algorithm of the RM.

## **4.2 Examples**

The working principle of the recursive method is illustrated by three examples. The setups for the first two academic examples only differ in the prescribed impact force. For an exponentially decreasing impact force, a solution of the inverse problem exists. Moreover, for this particular setup, an exact reference solution exists which is compared with the corresponding results obtained by applying the RM. With the second academic example, it is shown that it is only possible to calculate a solution up to a certain time if a linear increasing impact force is generated. The last example is devoted to the practical

application of percussion drilling. Thus, the impedance of the impacting piston is determined to generate a desired incident force in the rock.

### **4.2.1 Academic Examples**

The setup illustrated in Fig. 4.5 is considered where a semi-infinite rod (Rod 2) with constant impedance <sup>𝑍</sup>˜ is given. Rod 2 is at rest at the beginning (𝑣˜in <sup>=</sup> 0). Moreover, the velocity <sup>𝑣</sup>in of Rod 1 is known. The aim is to determine the unknown impedance function 𝑍(𝑠) of Rod 1 which generates a prescribed impact force 𝑁I(𝑡). The parameters

**Figure 4.5:** Setup of the academic examples. The aim is to determine the impedance 𝑍 which generates the impact force <sup>𝑁</sup>I(𝑡) for given velocity <sup>𝑣</sup>in and given impedance <sup>𝑍</sup>˜.



**Table 4.2:** Parameters for the academic examples.

#### **Example 1: Exponentially Decreasing Impact Force**

Initially, an exponentially decreasing impact force

$$N\_{\rm I}(t) = -\hat{N}\,\mathbf{e}^{\left(-\frac{t}{t\_0}\right)}\tag{4.24}$$

with arbitrary time constant <sup>𝑡</sup><sup>0</sup> is given. Lundberg and Lesser [47] have introduced the traveling time of the wavefront starting from the impact surface traveling to position 𝑠 with

$$\theta(\mathbf{s}) = \int\_0^s \frac{1}{c(\mathbf{s'})} \, \mathbf{ds'} \tag{4.25}$$

and derived the exact solution

$$Z\_{\text{exact}}(\theta) = \frac{Z\_0}{\left(1 + \frac{\tilde{Z} + Z\_0}{\tilde{Z}\iota\_0}\theta\right)^2}, \qquad Z\_0 = \frac{1}{\frac{v\_{\text{in}}}{\tilde{N}} - \frac{1}{\tilde{Z}}}.\tag{4.26}$$

Since the wavefront is traveling through one element within one time step, the impedances of the RM depending on the traveling time are expressed by

$$Z(\theta = q\Delta t) = Z\_q, \quad q = 0 \dots n-1. \tag{4.27}$$

The impedance results are depicted in Fig. 4.6. Both curves are continuously decreasing beginning from their maximum value at 𝜃(0). Physically, it is comprehensible that the

**Figure 4.6:** Comparison between the impedance function that is obtained with the RM and the exact solution.

curves decrease as the claimed impact force is decreasing exponentially as well. The impedance result determined with the RM matches very well the exact result. In fact, the maximum relative error is <sup>1</sup>.<sup>10</sup> · <sup>10</sup>−<sup>3</sup> %. The remaining difference to the exact solution is due to the piecewise constant approximation of the rod's impedances. The more elements are chosen the closer the results of the RM approach the exact solution.


An extract of the results at designated time steps is presented in Tab. 4.3. The simulation starts with the determination of the first element of Rod 1 and finishes with the last element (𝑞 <sup>=</sup> 100). From 𝑞 <sup>=</sup> <sup>10</sup> to until 𝑞 <sup>=</sup> 100, every tenth result is plotted.

**Table 4.3:** Simulation steps of the example with exponentially decreasing impact force.

There is a clear tendency to be observed: All values, that is the transmission factor <sup>𝐾</sup>𝑞 , the force difference <sup>Δ</sup>𝑁𝑞 , the force reserve 𝐾𝑞 𝑁<sup>0</sup>  and the impedances <sup>𝑍</sup>𝑞 are continuously decreasing. In Sec. 4.1.3, it has been proven that <sup>𝐾</sup>𝑞 is continuously decreasing. Therefore, the results are consistent with the theory. At the beginning of the simulation, the change in impedance is greatest. That is the reason, why <sup>𝐾</sup>𝑞 is decreasing strongest at that time. As a result of the decreasing factor <sup>𝐾</sup>𝑞 , the force reserves have to fall as well because <sup>𝑁</sup><sup>0</sup> is constant. The exponentially declining impact force leads to continuously declining impedances. The calculated values reveal that the presented example did make very little use of the force reserves.

#### **Example 2: Linearly Increasing Impact Force**

The second example only differs from the first one by the impact force which is now linearly increasing

$$N\_{\rm I}(t) = -\hat{N}\left(1 + 2\frac{t}{t\_0}\right). \tag{4.28}$$

In Fig. 4.7, the related impedance result is plotted on a logarithmic scale. At 𝑞 <sup>=</sup> 78, the simulation ends. Afterward, the calculated impedance value is negative which is

**Figure 4.7:** Impedance function of the example with linear increasing impact force.

not physical. Shortly before, the impedance values increase strongly. A more detailed overview of the results at designated simulation steps is presented in Tab. 4.4. The simulation starts with the first element of Rod 1 and finishes for 𝑞 <sup>=</sup> 78. From 𝑞 <sup>=</sup> <sup>10</sup> until 𝑞 <sup>=</sup> 70, every tenth result is plotted. At the time of interest, from 𝑞 <sup>=</sup> <sup>75</sup> on, every result is depicted. The absolute value of the force difference <sup>Δ</sup>𝑁𝑞  and the impedances <sup>𝑍</sup>𝑞 are continuously growing whereas the absolute value of the force reserve 𝐾𝑞 𝑁<sup>0</sup>  and the transmission factor <sup>𝐾</sup>𝑞 are continuously decreasing. The linear increasing impact


**Table 4.4:** Simulation steps of the example with linearly increasing impact force.

force leads to a continuous rise in the force difference and eventually to a continuous rise of the impedances. At the beginning, the increase is moderate, however, the longer the simulation takes the higher are the jumps in impedance. Especially, starting from 𝑞 <sup>=</sup> 70, sharp increases occur. For 𝑞 <sup>=</sup> 78, the existence condition is not met anymore. Therefore, no physical (positive) impedance can be determined.

One of the advantages of the presented approach is that the force condition immediately indicates if there will be a solution for the next time step. Moreover, it is possible to examine how long the given impact force can be generated. In this case, the maximum traveling time of the wavefront is <sup>77</sup>Δ𝑡.

### **4.2.2 Percussion Drilling**

Percussion drilling has a long tradition and is widely used for drilling applications. The classical setup (see Fig. 4.8) consists of a pneumatically or hydraulically actuated piston which hits onto the shank adapter at hitting frequencies between 50 and 100 Hz. The shank adapter together with the coupling sleeve transfer the kinetic energy of the piston into the rod. The resulting waves propagate through the drill rod to the drill bit and finally fracture the ground due to high impact forces. The emerging cuttings are flushed up by a high pressured fluid that is pumped to the ground through hollow rods. The rotation of the drilling rod at rotation speeds between 60 and 200 rpm leads to continuously changing points of impact [61]. Because of the energy losses during the transmission of the shock waves, percussion drilling is used for drilling lengths no longer than 30 m [37].

**Figure 4.8:** Sketch of the percussion drilling process, turned by 90◦ .

#### **Problem Description and Modeling**

The mechanics of percussion drilling has been studied for many years. However, the research has mainly focused on the direct problem. The modeling is either realized with a continuous model [9, 49], with an FEM model [13] or with lumped parameters [20, 42]. Comparisons between the results based on the 3D FEM and results based on the 1D theory show that the differences are very little so that the 1D theory is applicable [19, 48].

The considered inverse problem is depicted in Fig. 4.9. A piston of given length (0.5 m) and given initial velocity <sup>𝑣</sup>in is hitting onto the shank adapter. The shank adapter is modeled by a segment of constant impedance 𝑍˜ <sup>1</sup>, the coupling sleeve by a segment with impedance 𝑍˜ <sup>2</sup> and the section with constant impedance <sup>𝑍</sup>˜ <sup>3</sup> corresponds to the rod. The bit is composed of a constant impedance part 𝑍˜ <sup>4</sup> and a linear increase from <sup>𝑍</sup>˜ <sup>4</sup> to <sup>𝑍</sup>˜ <sup>5</sup> at the end. The bit-rock interaction during loading is modeled by a resistance force 𝐹 with coefficient 𝑘1. Its validity has been confirmed experimentally [17, 18, 39].

**Figure 4.9:** Example of a device for percussion drilling. The main parts are modeled as rod elements of different diameters. All dimensions are in m.

Lundberg and Collet [49] have shown that the efficiency of the drilling process is maximized if the incident force is of exponential shape

$$N\_{\rm I}(t) = -\hat{N}\,\mathrm{e}^{\rm C\prime I} \tag{4.29}$$

with an arbitrary initial force 𝑁<sup>ˆ</sup> and a constant 𝐶. The investigation motivates the realization of an exponential incident force by adjusting the impedance of the piston.

The simulation is run with the parameters listed in Tab. 4.5. The chosen coefficient <sup>𝑘</sup><sup>1</sup> and the constant 𝐶 are characteristic for a hard rock [49]. The simulation starts with


**Table 4.5:** Parameters used for the percussion drilling example.

the transformation of the incident force <sup>𝑁</sup><sup>i</sup> to the corresponding impact force <sup>𝑁</sup><sup>I</sup> . This enables applying the RM to calculate the piston's impedance function subsequently.

### **Transformation of the Incident Force to the Impact Force**

In Fig. 4.10a, the incident force described in Eq. (4.29) and the transformed impact force are depicted. For better comparability of the results, the incident force is time-shifted by the dead time <sup>𝑡</sup>dt which is the travel time from the impact surface to the bit-rock interaction.

The curve of the impact force is above the curve of the incident force. The reason is that the forces applied to the impact surface are amplified while traveling to the bit-rock interaction. The impact force curve exhibits three major jumps whose origins are depicted in Fig. 4.10b. The first jump occurs when the double reflected waves at the coupling sleeve reach the bit-rock interface. The second jump is significantly larger. It originates from the waves that are first reflected at the bit-rock interface, subsequently reflected at the transition between bit and rod and finally travel back to the bit-rock transition. These waves are compressive which is the reason why the magnitude of the impact force decreases afterward. The last jump is due to the double reflected waves between the impact surface and the coupling sleeve.

**Figure 4.10:** Transformation of the incident force to the impact force.

#### **Application of the Recursive Method**

As the incident force has now been transformed to the corresponding impact force, the RM is applied. In Fig. 4.11, the impedance over the distance 𝑠 (starting for 𝑠 <sup>=</sup> <sup>0</sup> at the hitting surface) of the optimized piston is shown. As expected, the impedance is increasing as the absolute value of the impact force is increasing as well. The three jumps also emerge in the impedance curve. This is sensible since the jumps in piston impedance lead to the desired jumps in impact force.

Using the RM, the results are determined exactly. This is confirmed by a comparison of the desired incident force with the incident force obtained by running the simulation with the determined impedance function of the piston. The maximum percentage error is <sup>1</sup>.<sup>13</sup> · <sup>10</sup>−<sup>12</sup> % which is in the order of machine precision. Moreover, it is possible to determine how long a solution of the inverse problem exists. For the chosen piston length, a solution exists. However, the longer the piston is, the higher are the

**Figure 4.11:** Calculated impedance function of the piston.

impedances at the end of the piston. At some point, no solution exists anymore. The huge difference between maximum and minimum impedance restricts the practical applicability. Besides, it is only possible to generate the desired incident force twice the travel time through the piston. Afterward, no elements are left to shape the curve. However, the percussion drilling example has shown that the RM is also applicable if the desired force is prescribed at a position other than the impact surface. Besides, the impedance function of the rod may vary and linear boundary conditions may be considered.

# **5 Experimental Investigation**

In this chapter, the 1D model is experimentally validated. After presenting the experimental setup with its components, preinvestigations are considered. Several aspects have to be examined in advance to guarantee reliable results. The most important factor is the repeatability of the experimental results as non-repeatable results have limited significance. The strains are measured by strain gauges (SGs) that are attached to different circumferential and axial positions of the rod. It is assessed if bending or the influence of the rod support results in different SG measurements. Subsequently, it is evaluated if the results can be scaled according to the impact velocity. The preinvestigations finish with the determination whether the piston bearings influence the results.

After the preinvestigations, two different experiments, which only differ from the cross sectional area variation of the piston, are conducted. Finally, the measurements are evaluated and compared with simulation results obtained by applying the wave propagation method.

## **5.1 Experimental Setup**

The test rig consists of a gun that pneumatically accelerates a piston which finally hits the rod (see Fig. 5.1). The impact velocity of the piston is measured by a photoelectronic fork sensor called *SpeedGate* of the manufacturer *Frederikson*. The sensor's two vertical laser rays at known axial distance are cut by the front face of the piston at different times. As both the time difference of these cuts and the related axial distance is known, the impact velocity can be determined. After the hit, force waves propagate into piston and rod. At three different axial positions, SGs, that are only able to detect longitudinal elongation, are attached to the rod. Their gauge length (0.<sup>2</sup> mm) is very short which ensures good spatial resolution.

#### 5 Experimental Investigation

**Figure 5.1:** (a) Experimental setup and (b) its schematic sketch.

Each SG is part of a quarter bridge (see Fig. 5.2) that transforms the mechanical strain into electrical voltage. The transformed electrical signal is amplified and filtered by an analog anti-aliasing filter (Bessel, 5th order) with corner frequency 100 kHz. Subsequently, the signal is digitalized and filtered by a digital IRR filter. The filter properties of the analog filter are predefined whereas the properties of the digital filter are set according to the application. For the given example, a conflict of interest exists. On the one hand, the corner frequency of the digital filter (Bessel filter 4th order) has to be very high as the

**Figure 5.2:** Measuring chain from the SG to the representation at the computer.

transient strain signal includes high frequencies. On the other hand, the higher the corner frequency that is chosen the more noise is included in the signal. Based on these conflicting goals, the corner frequency 40 kHz is set. The signal processing is conducted with the signal conditioner *Sirius HS* from *ZSE*. Finally, the measurements are evaluated by a computer.

The high-speed camera records the hit at its maximum frame rate (38565 fps) which enables detecting at what time piston and rod separate. Moreover, the videos provide very good insight during the impact. As the ambient light is not sufficient for these frame rates, external LED-lights illuminate the impact zone. The parameters of the SGs, the high speed camera and the signal conditioner are summarized in Tab. 5.1.

The rod can move freely at its right end. The stop only prevents the rod from crashing into the fence after measuring. A more detailed examination of the piston and rod and of the pressure chamber is provided in the next subsections.


**Table 5.1:** Parameters of the SGs, the high speed camera and the signal conditioner.

#### **Piston and Rod**

For the experimental results, one rod and two different pistons are applied. The dimensions of the pistons are shown in Fig. 5.3 and Fig. 5.4, respectively. Piston 1 is symmetrical with two small diameter parts (40 mm) at their ends and one large diameter part (50 mm) in the middle. Piston 2 consists of many cross section area parts, including several steps and one cone.

The dimensions of the rod and the axial strain gauge positions (SGPs) are depicted in Fig. 5.5. The rod's diameter is constant (50 mm) with a total length of 1800 mm. The axial SGP that is closest to the impact zone is named SGP1, the one in the middle of the rod SGP<sup>2</sup> and the one to the very end of the rod is named SGP3. The SGs are attached symmetrically meaning that the distance of SGP<sup>1</sup> from the left end equals the distance of SGP<sup>3</sup> from the right end (each 110 mm). Furthermore, SGP<sup>2</sup> is exactly in the middle of the rod. The SGs are mounted at different circumferential positions. They are either

**Figure 5.3:** Dimensions of Piston 1 in mm.

**Figure 5.4:** Dimensions of Piston 2 in mm.

mounted on the top (0 ◦ ), right (90◦ ), bottom (180◦ ) or left (270◦ ) side of the rod (see Fig. 5.5). Overall, eight SGs are attached to the rod. Two SGs are attached to SGP<sup>1</sup> (top, bottom), four SGs are attached to SGP<sup>2</sup> (top, right, bottom, left) and two SGs are attached to SGP<sup>3</sup> (top, bottom). The SGs are named by SGi,x, where i refers to the SGP and x to the angle. For example, SG1,90◦ refers to the right SG at SGP1.

**Figure 5.5:** Dimensions of the rod and the axial position of the SGs in mm (left) and a sketch of the circumferential positions of the SGs (right).

The pistons and the rod are made of steel with the parameters listed in Tab. 5.2.


**Table 5.2:** Parameters for the pistons and the rod.

#### **Pressure Chamber**

The acceleration of the piston is carried out pneumatically. Between the outer and inner diameter of the gun, a pressure chamber is filled with pressurized air (see Fig. 5.6). The chamber is kept closed by a bushing which is actuated pneumatically until the pressure in the chamber reaches the desired air pressure. Once the pneumatic actuator pulls the

**Figure 5.6:** Pressure chamber of the gun in (a) closed and (b) open state.

bushing to the left, the pressurized air is flowing behind the left end of the piston which leads to its acceleration. The experimental setup is a single hit device meaning that the piston is not automatically moving back and forth. Moreover, as the pressure in the pressure chamber is controlled, it is not possible to control the piston's impact velocity directly.

## **5.2 Experimental Procedure**

Two different setups are investigated. The first setup (see Fig. 5.7) only differs from the second setup (see Fig. 5.8) in the shape of the impacting piston. The rod, which initially is at rest, can move freely at its right end. The stop does not influence the results as only the transient behavior of the waves immediately after the impact is analyzed. The initial velocity of the pistons may be changed according to the investigated example.

**Figure 5.7:** Setup 1.

**Figure 5.8:** Setup 2.

For the evaluation of the results, the normal forces are calculated from the measurements at the SGs. The normal force <sup>𝑁</sup>𝑖,𝑥 refers to the force at SGi,x. For example, <sup>𝑁</sup>2,<sup>0</sup> ◦ is the normal force at SG2,0 ◦ . For most evaluations, only the normal force <sup>𝑁</sup>𝑘 , which compensates bending, is denoted. It is determined with

$$N\_k = \frac{1}{2} \left( N\_{k,0^\circ} + N\_{k,180^\circ} \right), \ k = 1 \ldots 3. \tag{5.1}$$

The subscript number of the forces refers to the corresponding SGP.

The SGs are attached to the top and bottom at each SGP to both compensate the normal forces for bending and to determine the maximum bending stress

$$
\sigma\_{k, \text{bend}} = \frac{1}{2 \, A\_{\text{Rod}}} \left( N\_{k, 0^\circ} - N\_{k, 180^\circ} \right), \text{ } k = 1 \, \dots \text{3}, \tag{5.2}
$$

with the cross section area <sup>𝐴</sup>Rod of the rod. By using two diametrically opposite strain gauges in a half bridge, it would not be possible to determine the maximum bending stresses which is why a quarter bridge has been used.

In Tab. 5.3, an overview of the conducted experiments is given. The preinvestigations are all carried out with the first setup. The *repeatability* experiment is conducted three times. Since it is not possible to directly control the impact velocity, little differences in



**Table 5.3:** Parameters of the conducted experiments.

velocity occur. However, the deviations are negligible. The corresponding results for the normal force at SGP<sup>1</sup> are compared. The *circumferential SGP* experiments aim at both determining the influence of the bending stress and the influence of the rod bearings. Thus, the four normal forces at SGP<sup>2</sup> are compared and the ratio max <sup>𝜎</sup>2,bend /max <sup>|</sup>𝜎2<sup>|</sup> with <sup>𝜎</sup><sup>2</sup> <sup>=</sup> <sup>𝑁</sup>2/𝐴Rod is determined. The *impact velocity* experiments differ in the impact velocity which range from <sup>3</sup>.<sup>47</sup> m s−<sup>1</sup> to <sup>5</sup>.<sup>60</sup> m s−<sup>1</sup> . It is assessed if the results are scalable by comparing the normalized normal forces at SGP1. The *piston bearings* experiments give information on whether the piston bearings have to be considered in the model by comparing the results for 𝑁1.

In the main investigations section, the normal forces <sup>𝑁</sup>1, 𝑁<sup>2</sup> and <sup>𝑁</sup><sup>3</sup> are analyzed and interpreted both for Setup 1 and Setup 2. At the end, the experimental results for <sup>𝑁</sup><sup>1</sup> are compared with the corresponding simulation results which are obtained by applying the wave propagation method.

## **5.3 Preinvestigations**

The focus of the preinvestigations is on determining possible influence factors on the measurements. Therefore, the interpretation of the results is qualitatively and comparative. A detailed interpretation of the results will be conducted in Sec. 5.4.

### **5.3.1 Repeatability**

In Fig. 5.9, the experimental results at SG1,0 ◦ for three impacts are depicted. The results

**Figure 5.9:** Examination of the repeatability of the mesaurement results. The initial velocities of the piston are in the range between <sup>4</sup>.57 m s−<sup>1</sup> and <sup>4</sup>.60 m s−<sup>1</sup> .

match very well without any obvious deviations. Very similar results are achieved when the measurements are analyzed at other axial and other circumferential positions. Therefore, it is concluded that the results are repeatable.

### **5.3.2 Influence of Circumferential Strain Gauge Position**

The SGs are placed at different circumferential positions. In theory, a direct central impact only causes longitudinal waves to propagate into the rod. However, optimal conditions are never met. Furthermore, the rod is supported at two positions which may influence the results due to friction. This examination aims to assess if the rod bearings have to be considered in modeling and if the impact is direct central. Therefore, the results at the second SG for four different circumferential positions are shown in Fig. 5.10.

**Figure 5.10:** Influence of the circumferential SGP on the experimental results. The initial velocity of the piston is <sup>4</sup>.60 m s−<sup>1</sup> .

The qualitative behavior of the results is almost the same for all angles. Quantitatively, the deviations are very small which can be regarded in the zoom window. Therefore, it is concluded that the influence of the rod bearing is negligible and that the impact is close to direct central. This assumption is confirmed by the calculated ratio

$$e = \frac{\max \left| \sigma\_{2,\text{bend}} \right|}{\max |\sigma\_2|} 100 \tag{5.3}$$

between maximum bending stress and normal stress which is <sup>2</sup>.74 %.

### **5.3.3 Influence of Impact Velocity**

The initial velocity of the impacting piston varies for different applications. As the presented test rig is also restricted in terms of its maximum piston velocity, the question arises if the results are scalable. In Fig. 5.11, the measurement results obtained with three different initial velocities are depicted. The results are normalized by the absolute value of the respective global minimum. The deviations between the normalized results are very low which is highlighted in the zoom box. This allows the conclusion that within the presented velocity range from <sup>3</sup>.<sup>47</sup> m s−<sup>1</sup> until <sup>5</sup>.<sup>60</sup> m s−<sup>1</sup> the results are scalable. However, it is expected that the results are also scalable for a wider range but this needs finally to be proven by other experiments.

**Figure 5.11:** Influence of the initial velocity of the piston on the experimental results.

### **5.3.4 Influence of Piston Bearings**

The pistons are mounted with bearings made of Polytetrafluorethylen as can be seen in Fig. 5.12. Since it is not possible to accelerate the piston without bearings with the presented test rig configuration, it has to be clarified if the piston bearings influence the measurements. If the bearings influenced the results, there would be a difference in the measurements depending on the bearing width. However, as can be seen in Fig. 5.12,

**Figure 5.12:** Influence of the bearings of the piston on the experimental results. The initial velocities of the piston are in the range between <sup>4</sup>.60 m s−<sup>1</sup> and <sup>4</sup>.62 m s−<sup>1</sup> .

there is no significant difference apparent which is why the bearings are not considered.

In summary, the preinvestigations reveal that the experimental results are


## **5.4 Main Investigations**

The goal of the main investigations is to validate the 1D model presented in Chap. 3. Therefore, two experiments are conducted, analyzed and finally compared with the related simulation results.

The first experiment is conducted to assess if the experimental results are reliable. As the geometry of Piston 1 varies only very little, it is possible to compare the measurements with analytical results.

The second experiment is conducted with a complex piston geometry (Piston 2) with many changes of the cross section area in axial direction. Piston 2 is applied to evaluate if the 1D model is also valid for complex piston shapes.

### **5.4.1 Standard Impacting Rod**

The first experiment coincides with the setup used for the preinvestigations (see Fig. 5.7). Piston 1 is hitting the rod at initial velocity <sup>4</sup>.60 m s−<sup>1</sup> .

### **Experimental Results**

The time development of the forces, that are related to the corresponding SGs at three different SGPs, is illustrated in Fig. 5.13. The numbers on the graphs belong to selected time points. The qualitative behavior at these time points is explained with Fig. 5.14. At 𝑡1, the first deflection is observed at the first SGP. After an abrupt drop of the force, the force remains almost constant for a short period. Subsequently, another step to the maximum compressive force follows. Apart from small oscillations after 𝑡3, the wave shape between <sup>𝑡</sup><sup>1</sup> and <sup>𝑡</sup><sup>3</sup> is symmetrical. Just before <sup>𝑡</sup>1, the piston hits the left end of

**Figure 5.13:** Forces versus time (first experiment) gained at the SGPs. The initial velocity of the piston is <sup>4</sup>.60 m s−<sup>1</sup> .

the rod. The time from <sup>𝑡</sup> <sup>=</sup> <sup>0</sup> until <sup>𝑡</sup><sup>1</sup> arises due to the pretrigger settings and shall not be confused with the traveling time of the waves from the hitting surface to the first SGP. During the time from <sup>𝑡</sup><sup>3</sup> to <sup>𝑡</sup>8, the force is close to zero. Subsequently, a point symmetrical section with a tensile force peak at the beginning and a compressive force peak at the end (between <sup>𝑡</sup><sup>8</sup> and <sup>𝑡</sup>10) is observed.

The incident wave shape at the second SGP from <sup>𝑡</sup><sup>2</sup> to <sup>𝑡</sup><sup>5</sup> is very similar but time-delayed with respect to the wave shape at the first SGP between <sup>𝑡</sup><sup>1</sup> and <sup>𝑡</sup>3. Furthermore, the shape of the subsequent part (𝑡<sup>6</sup> <sup>−</sup> <sup>𝑡</sup>9) resembles the incident wave with the opposite sign. The sequence from <sup>𝑡</sup><sup>2</sup> to <sup>𝑡</sup><sup>9</sup> also repeats shortly after.

The results at the third SGP lastly drop (starting from 𝑡4). The following wave shape between <sup>𝑡</sup><sup>4</sup> and <sup>𝑡</sup><sup>7</sup> resembles the wave shape of the first SGP between <sup>𝑡</sup><sup>8</sup> and <sup>𝑡</sup>10. However, this time the wave shape begins with a compressive force peak and ends with a tensile peak which is reverse at the first SGP. Again, the wave shape between <sup>𝑡</sup><sup>4</sup> and <sup>𝑡</sup><sup>7</sup> repeats some time later.

For understanding the qualitative behavior of the force waves, it is advantageous to study the hit of two rods of the same impedance. In Fig. 5.14, the forces at any piston and rod position are depicted for selected time points which refer to the time points marked in Fig. 5.13. The piston is located on the left side whereas the rod with its three SGPs is placed on the right side. The force wave propagation direction is indicated with small arrows and the colors of the waves state if a wave is compressional (green) or tensile (orange).

After the hit, compressional force waves travel both into piston and rod. At 𝑡1, the wavefront reaches the first SGP and leads to a jump. It is impossible to perfectly represent a jump in the SG measurements due to four reasons. Firstly, the SGs always measure the average elongation along their gauge length. Secondly, unlimited sample rates cannot be chosen. Thirdly, the signal is filtered with a low pass filter. And most importantly, due to 3D-effects, there is always a finite rise in reality. Moreover, the filter causes oscillations after jumps. This has to be considered when interpreting the results.

At 𝑡2, the wavefront of the rod reaches the second SGP leading to the same jump as at the first SGP. Meanwhile, the wavefront of the piston has traveled to the left end where the wave is reflected as a tensile force. The compressive and tensile force waves are equal in amount which is why the resulting force is zero. The tensile force wave continues

5.4 Main Investigations

**Figure 5.14:** Schematic sketch of the wave propagation after a hit of two rods with same impedance.

traveling along the piston until it reaches the transition to the rod. This is the time <sup>𝑡</sup><sup>s</sup> when piston and rod separate as no tensile force may be transmitted at the transition. Therefore, a rectangular wave shape that extends to double the length of the piston, propagates along the rod. The development of the forces in the piston is not pursued any longer.

At 𝑡3, the compression wave passes the first SGP leading to zero force again. The time period 𝑡 <sup>∗</sup> between the arrival of the wave front at the first SGP until the wave has completely passed the first SGP is determined with

$$t^\* = \frac{2\ell\_{\rm Pis}}{c}.\tag{5.4}$$

For the given example, 𝑡 ∗ equals <sup>1</sup>.<sup>93</sup> · <sup>10</sup>−<sup>3</sup> s. This is in very good accordance with the time span between <sup>𝑡</sup><sup>3</sup> and <sup>𝑡</sup>1.

As the rod's impedance is constant, the wave shape does not change while traveling along the rod. Therefore, the wavefront, that arrives at the third SGP at time 𝑡4, has the same shape when arriving at the other SGPs. At SGP2, no reflecting waves interfere which is why the shape of the normal force at the second SGP (between <sup>𝑡</sup><sup>2</sup> and <sup>𝑡</sup>5) is the time-delayed shape of the force at the first SGP (between <sup>𝑡</sup><sup>1</sup> and <sup>𝑡</sup>3).

At 𝑡5, two events happen almost simultaneously. The incident compressive wave passes the second SGP and the reflected tensile force wavefront reaches the third SGP leading at both SGPs to zero resulting force.

When the compressive wave passes the last SGP at 𝑡6, only the tensile force is acting on the SGs which is the reason why the force is jumping again. At the same time, the tensile wavefront reaches the second SGP which causes a sharp rise in tensile force.

When the reflected part of the incident wave completely passes the third SGP (time 𝑡7), the force drops to zero. From time <sup>𝑡</sup><sup>8</sup> to <sup>𝑡</sup>10, the same procedure as for the third SGP takes place for the first SGP with the only difference that the incident wave is tensile initially.

At <sup>𝑡</sup>10, the state equals the state at <sup>𝑡</sup>3. Thus, the period between <sup>𝑡</sup><sup>3</sup> and <sup>𝑡</sup><sup>10</sup> is repeated all over again until the oscillations are damped by, e.g., material damping. However, in the considered period, the influence of material damping is negligible. The period time complies with the time it takes to travel twice along the rod.

#### **Comparison with Numerical Results**

The measurements are compared with the simulation results which are generated by applying the wave propagation method. The parameters, the initial conditions and the setup of the simulation are set identically to the experiment. Moreover, the results are filtered by the same filter (Bessel 4th order) as the measurements are filtered. The

**Figure 5.15:** Evaluation of the first experiment.

resulting graph is plotted in Fig. 5.15a. A part of the first incident wave is highlighted in the zoom box. The deviations between the results are very little. The largest deviations are detected after jumps when the experimental solution still oscillates. Both the qualitative and the quantitative behavior of the curves matches very well as can be seen in detail in the zoom box. After the arrival of the incident compressive force, four main deflections are observed. Its sources are schematically visualized in Fig. 5.15b.

The first deflection originates from the wave that incidents in the piston, reflects at its first step and finally transmits to SGP1. The second deflection is considerably smaller which is due to more reflections included in the development. The waves responsible for the first and second step are identical until the wavefront reaches the transition between piston and rod. While the first wave is transmitted to the SGs, the second wave is reflected back and forth twice between transition and the first step of the piston before it arrives at the SGs. The third wave is transmitted to the second step of the piston where it is reflected and transmitted directly to the SGs. The fourth wave is similar to the third one with the difference that the fourth wave is reflected at the free end of the piston. It is important to emphasize that the presented waves are only the ones that contribute most to the final shape of the force wave. The other waves are reflected very often so that they are very low in amplitude and negligible for the qualitative understanding of the final wave shape.

In summary, the accordance of the simulation with the experimental results is very good. Even small deflections are depicted by both approaches. Similar results are achieved at the remaining two SGPs. Therefore, the 1D model is validated for small variations in the piston geometry. In the next subsection, the question is answered if the statement is also valid for major variations of the geometry.

### **5.4.2 Complex Impacting Rod**

The second experiment is similar to the first experiment but with a different piston hitting the rod (see Fig. 5.8). Piston 2 has several cross sectional area variations, including a conical increase at the beginning.

### **Experimental Results**

In Fig. 5.16, the experimental results of the second experiment at three SGPs are presented. Since only the piston geometry varies compared to the first experiment, the qualitative behavior of the results is similar. An incident wave starting at 𝑡 <sup>≈</sup> <sup>1</sup> ms at the first SGP is traveling along the rod maintaining its shape while passing the second SGP. Over time, the incident wave shape continuously alternates between compressive and tensile force at SGP2. At the third SGP, only the very beginning of the incident wave is visible as the measurement is superimposed by reflecting waves at the rod's free end a little while later. After two peaks (compressive and tensile), a section with forces close to zero follows before the signal is repeated. Very similar results occur at the first SGP after

**Figure 5.16:** Forces versus time (second experiment) gained at the SGPs. The initial velocity of the piston is <sup>4</sup>.58 m s−<sup>1</sup> .

the incident wave has passed with the difference that the sequence starts with a tensile peak.

### **Comparison with Numerical Results**

The comparison of the experimental versus the simulation results at the first SGP is presented in Fig. 5.17a. For detailed analysis, the incident wave is magnified in the zoom box (filtered on the left side and unfiltered on the right side of Fig. 5.17b). It is only possible to generate filtered results for the simulation as unfiltered measurements would include aliasing effects. The final shape of the incident wave consists of a multitude of superimposed waves that are transmitted and reflected several times within the piston. However, some waves significantly contribute to the outcome. These waves are illustrated in Fig. 5.17c.

The qualitative and quantitative behavior matches very well with only small deviations between the results. The greatest qualitative deviation occurs approximately in the middle of the filtered results. At this position, the simulation result is shortly deflecting in negative direction while the experimental result is continuously increasing. Rod and piston do not separate when the reflected wave at the free end of the piston reaches the transition zone but shortly later which leads to a wider incident wave compared to the first experiment. It is very difficult to analyze the origin of every oscillation in the results as the longer the time proceeds the more piston variations influence the force shape. Moreover, the steps blur due to filtering. That is the reason why unfiltered results are used to analyze the most significant steps.

The first step at <sup>𝑡</sup><sup>1</sup> originates from the wave which propagates into the piston, reflects at the beginning of the cone and transmits back to the first SGP. From <sup>𝑡</sup><sup>1</sup> until close to 𝑡2, the force decreases almost linearly which is due to the conical part of the piston. When the wave, which reflects at the end of the cone, reaches SGP<sup>1</sup> at <sup>𝑡</sup>2, the force rises immediately. The steps at <sup>𝑡</sup>3/𝑡<sup>4</sup> refer to the waves that are reflected at the beginning/end of the middle block. The reflected wave at the small block close to the end of the piston is responsible for the jump in force at 𝑡5. The biggest jump in force occurs at 𝑡6. It originates from the wave that transmits through the piston, reflects at the piston's free end and transmits back to SGP1.

While comparing the unfiltered with the filtered results, it is noticeable that not every detail is represented in the filtered result. However, the main changes are visible. The comparison shows that the main variations in the piston geometry are depicted in both

**Figure 5.17:** Evaluation of the second experiment.

experimental and simulation results with very good accordance. Therefore, it can be concluded that the 1D model is validated. Moreover, the comparisons suggest that for the given impact problems, the 1D consideration is sufficient. The advantage of the 1D theory is that the influence of the respective piston sections may be worked out. Furthermore, the calculation time is much shorter. The main advantage is, however, that it is possible to develop a method that solves the inverse problem based on the 1D theory.

# **6 Conclusions**

The thesis aims to develop a method related to inverse problems in impacting rods. Specifically, the impedance of the impacting rod is determined so that a prescribed force at an arbitrary axial position of the second rod is generated (*impactor synthesis problem*). The method is relevant for several experimental and practical applications. For example, in percussion drilling, the question arises how the impedance of the impacting rod can be determined to generate an incident force that improves the efficiency of the drilling process.

The *impactor synthesis problem* is only addressed by a few researchers. The most promising approach solves the inverse problem in terms of an integral equation. By applying the integral based approach, several approximation methods are carried out to calculate the impactor's impedance. Moreover, it is only possible to prescribe the desired force at the hitting surface. In the context of the thesis, an approach based on the idea of differential methods is developed. It addresses the drawbacks of the integral based method mentioned above.

The method relies on the idea of directly exploiting the structure of the wave propagation problem. Therefore, the impedance functions of the rods are approximated by piecewise constant step functions. Use is made of the circumstance that each element of the unknown impactor is influencing the prescribed force at different moments for the first time. The element, which just influences the actual force for the first time, is adjusted so that the actual force matches the prescribed force. Thus, the elements are determined recursively starting from the element which is closest to the prescribed force position.

For piecewise constant impedance functions, the developed method delivers exact solutions in closed-form. Moreover, a condition is derived which states if a solution of the inverse problem exists. It is shown that any additional element, whose impedance differs from the previous one, diminishes the probability that a solution exists. The method allows to prescribe the force at any axial position of the second rod. Furthermore, arbitrary initial conditions are applicable to the second rod, including external forces.

The practical applicability of the method is validated in two steps. Firstly, the determined impedance function generates the prescribed force over time which is calculated by applying the 1D model. Secondly, the 1D model is validated experimentally by comparing experimental versus simulation results. In particular, it is shown that the model is also valid for impacting rods with complex geometry variations. As scattering problems are special cases of the *impactor synthesis problem*, the developed method is widely applicable.

# **List of Abbreviations**


# **List of Symbols**

#### **Indices**


### **Symbols**




# **Bibliography**


# **Publications**


#### ISSN 1614-3914


### Band 2 Clemens Reitze Closed Loop, Entwicklungsplattform für mechatronische Fahrdynamikregelsysteme. 2004 ISBN 3-937300-19-8

### Band 3 Martin Georg Cichon Zum Einfluß stochastischer Anregungen auf mechanische Systeme. 2006 ISBN 3-86644-003-0

### Band 4 Rainer Keppler Zur Modellierung und Simulation von Mehrkörpersystemen unter Berücksichtigung von Greifkontakt bei Robotern. 2007 ISBN 978-3-86644-092-0

### Band 5 Bernd Waltersberger Strukturdynamik mit ein- und zweiseitigen Bindungen aufgrund reibungsbehafteter Kontakte. 2007 ISBN 978-3-86644-153-8

### Band 6 Rüdiger Benz

Fahrzeugsimulation zur Zuverlässigkeitsabsicherung von karosseriefesten Kfz-Komponenten. 2008 ISBN 978-3-86644-197-2

### Band 7 Pierre Barthels

Zur Modellierung, dynamischen Simulation und Schwingungsunterdrückung bei nichtglatten, zeitvarianten Balkensystemen. 2008 ISBN 978-3-86644-217-7

### Band 8 Hartmut Hetzler

Zur Stabilität von Systemen bewegter Kontinua mit Reibkontakten am Beispiel des Bremsenquietschens. 2008 ISBN 978-3-86644-229-0

Band 9 Frank Dienerowitz Der Helixaktor – Zum Konzept eines vorverwundenen Biegeaktors. 2008 ISBN 978-3-86644-232-0

### Band 10 Christian Rudolf Piezoelektrische Self-sensing-Aktoren zur Korrektur

statischer Verlagerungen. 2008 ISBN 978-3-86644-267-2

### Band 11 Günther Stelzner

Zur Modellierung und Simulation biomechanischer Mehrkörpersysteme. 2009 ISBN 978-3-86644-340-2

### Band 12 Christian Wetzel

Zur probabilistischen Betrachtung von Schienen- und Kraftfahrzeugsystemen unter zufälliger Windanregung. 2010 ISBN 978-3-86644-444-7

### Band 13 Wolfgang Stamm

Modellierung und Simulation von Mehrkörpersystemen mit flächigen Reibkontakten. 2011 ISBN 978-3-86644-605-2

### Band 14 Felix Fritz

Modellierung von Wälzlagern als generische Maschinenelemente einer Mehrkörpersimulation. 2011 ISBN 978-3-86644-667-0

### Band 15 Aydin Boyaci

Zum Stabilitäts- und Bifurkationsverhalten hochtouriger Rotoren in Gleitlagern. 2012 ISBN 978-3-86644-780-6

### Band 16 Rugerri Toni Liong Application of the cohesive zone model to the analysis of rotors with a transverse crack. 2012 ISBN 978-3-86644-791-2

### Band 17 Ulrich Bittner

Strukturakustische Optimierung von Axialkolbeneinheiten. Modellbildung, Validierung und Topologieoptimierung. 2013 ISBN 978-3-86644-938-1

### Band 18 Alexander Karmazin

 Time-efficient Simulation of Surface-excited Guided Lamb Wave Propagation in Composites. 2013 ISBN 978-3-86644-935-0

### Band 19 Heike Vogt

Zum Einfluss von Fahrzeug- und Straßenparametern auf die Ausbildung von Straßenunebenheiten. 2013 ISBN 978-3-7315-0023-0

### Band 20 Laurent Ineichen

Konzeptvergleich zur Bekämpfung der Torsionsschwingungen im Antriebsstrang eines Kraftfahrzeugs. 2013 ISBN 978-3-7315-0030-8

### Band 21 Sietze van Buuren

Modeling and simulation of porous journal bearings in multibody systems. 2013 ISBN 978-3-7315-0084-1

### Band 22 Dominik Kern

Neuartige Drehgelenke für reibungsarme Mechanismen. 2013 ISBN 978-3-7315-0103-9

### Band 23 Nicole Gaus

Zur Ermittlung eines stochastischen Reibwerts und dessen Einfluss auf reibungserregte Schwingungen. 2013 ISBN 978-3-7315-0118-3

### Band 24 Fabian Bauer

 Optimierung der Energieeffizienz zweibeiniger Roboter durch elastische Kopplungen. 2014 ISBN 978-3-7315-0256-2

### Band 25 Benedikt Wiegert

Nichtlineare Schwingungen von Systemen mit elastohydrodynamischen Linienkontakten. 2015 ISBN 978-3-7315-0350-7

### Band 26 Arsenty Tikhomolov

Analytische, numerische und messtechnische Untersuchung der Dynamik von Fahrzeugkupplungen am Beispiel des Trennproblems. 2015 ISBN 978-3-7315-0362-0

### Band 27 Daniel Maier

 On the Use of Model Order Reduction Techniques for the Elastohydrodynamic Contact Problem. 2015 ISBN 978-3-7315-0369-9

### Band 28 Xiaoyu Zhang

Crosswind stability of vehicles under nonstationary wind excitation. 2015 ISBN 978-3-7315-0376-7

### Band 29 Jens Deppler

Ein Beitrag zur viskoelastischen Modellierung nichtholonomer Bindungsgleichungen. 2017 ISBN 978-3-7315-0548-8

### Band 30 Georg Jehle

Zur Modellbildung und Simulation reibungserregter Schwingungen in Pkw-Schaltgetrieben. 2017 ISBN 978-3-7315-0668-3

### Band 31 Joachim Klima

Lubricant transport towards tribocontact in capillary surface structures. 2018 ISBN 978-3-7315-0814-4

### Band 32 Gábor Kenderi

Nonparametric identification of nonlinear dynamic systems. 2018 ISBN 978-3-7315-0834-2

#### Band 33 Jan Henrik Schmidt

An efficient solution procedure for elastohydrodynamic contact problems considering structural dynamics. 2018 ISBN 978-3-7315-0872-4

### Band 34 Ulrich Johannes Römer

Über den Einfluss der Fußgeometrie auf die Energieeffizienz beim zweibeinigen Gehen. 2019 ISBN 978-3-7315-0887-8

### Band 35 Simon Kapelke

Zur Beeinflussung reibungsbehafteter Systeme mithilfe überlagerter Schwingungen. 2019 ISBN 978-3-7315-0905-9

### Band 36 Jens Burgert

On direct and inverse problems related to longitudinal impact of non-uniform elastic rods. 2021 ISBN 978-3-7315-1087-1

Impacting rods are used in many practical and experimental applications. The analysis problem of finding the time-dependent impact force for given properties of the impacting rods and known initial velocities is thoroughly studied in literature. However, very little literature has been published which addresses its inverse synthesis problem: Find the location-dependent impedance function of the impacting rod that generates a prescribed impact force.

In this contribution, a method that solves the synthesis problem is developed. The impedances of the linear elastic rods are discretized by elements of piecewise constant impedance. The method utilizes that each element of the unknown impactor influences the prescribed force at different time instants for the first time. Thus, the unknown impedances of the impactor may be determined recursively. The developed method delivers exact solutions in closed-form. Moreover, a condition is derived which states if a physically meaningful solution exists.

Finally, the underlying 1D model is validated experimentally. To this end, a single-hit test rig is set up. The comparison of simulation versus experimental results yields good accordance both for standard and complex geometries of the impacting rods.

On direct and inverse problems related to longitudinal impact of non-uniform elastic rods

Jens Burgert

ISSN 1614-3914 ISBN 978-3-7315-1087-1

Gedruckt auf FSC-zertifiziertem Papier